cannin/biostars_qa huggingface dataset:
Question: How to Grep the complete sequences containing a specific motif in a fasta file using shell command? Also, I want to include the lines beginning with a `>` before these target sequences.

I found this post :https://www.biostars.org/p/274859/ similar to my problem but I'm looking for different motif:

My motifs look like that :

SXXXX(F/S)XXXL

All my fasta file in one line and I have more than 300 sequences:

for example:

my sequence :

>sp|Q9H257.2|CARD9_HUMAN RecName: Full=Caspase recruitment domain-containing protein 9; Short=hCARD9
MSDYENDDECWSVLEGSRVTLTSVIDRSRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9H37.2|CTYU_HUMAN
HHHSVLEGFRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9re7.2|CARer_HUMAN RecName
BKLSVLEGWRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

The result should be only the first two sequences because they have the motifs `SXXXX(F/S)XXXL`

>sp|Q9H257.2|CARD9_HUMAN RecName: Full=Caspase recruitment domain-containing protein 9; Short=hCARD9
MSDYENDDECWSVLEGSRVTLTSVIDRSRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9H37.2|CTYU_HUMAN
HHHSVLEGFRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

I tried these command but that returned all three sequences

grep 'S...F\|S\|L.\(.\)\1\{4\}' jara3.fasta -B 1 > jara4.fasta

Answer: using [seqkit][1], for sequences with regex (input from OP):

$ seqkit grep -srip 'S.{4}[FS].{3}L' test.fa

>sp|Q9H257.2|CARD9_HUMAN RecName: Full=Caspase recruitment domain-containing protein 9; Short=hCARD9
MSDYENDDECWSVLEGSRVTLTSVIDRSRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRK
VGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLM
TEVMKLQKKVQDLTALLSSK
>sp|Q9H37.2|CTYU_HUMAN
HHHSVLEGFRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDIL
QRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQK
KVQDLTALLSSK

For sequences without regex:

$ seqkit grep -svrip 'S.{4}[FS].{3}L' test.fa

>sp|Q9re7.2|CARer_HUMAN RecName
BKLSVLEGWRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDIL
QRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQK
KVQDLTALLSSK

Note: I used . in sequence assuming that only alphabets (AA) are present in sequence.
[1]: https://bioinf.shenwei.me/seqkit/
Question: HI

I am trying to parse a fasta file by a key word in the description using biopython (1.73). I am a novice and perhaps this question has been asked and answered but I cant seem to find the solution myself.

My fasta file looks like this. I want to create a file that contains only the sequence with gene=balVIM

>lcl|KP771862.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]
ATGTTCAAACTTTTGAGTAAGTTATTGGTCTATTTGACCGCGTCTATCATGGCTATTGCGAGTCCGCTCG
CTTTTTCCGTAGATTCTAGCGGTGAGTATCCGACAGTCAGCGAAATTCCGGTCGGGGAGGTCCGGCTTTA
>lcl|KP749829.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]
ATGTTCAAATTTTTGAGTAAGTTATTGGTCTATTTGACCGCGTCTATCATGGCTATTGCGAGTCCGCTCG
CTTT
>lcl|EU118148.2_gene_1 [gene=intI1] [location=complement(<1..231)] [gbkey=Gene]
ATGAAAACCGCCACTGCGCCGTTACCACCGCTGCGTTCGGTCAAGGTTCTGGACCAGTTGCGTGAGCGCA
TACGCTACTTGCATTACAGTTTACGAACCGAACAGGCTTATGTCAACTGGGTTCGTGCCTTCATCCGTTT
CCACGGTGTGCGTCACCCGGCAACCTTGGGCAGCAGCGAAGTCGAGGCATTTCTGTCCTGGCTGGCGAAC
GAGCGCAAGGTTTCGGTCTCC
>lcl|EU118148.2_gene_2 [gene=aacA29a] [location=485..880] [gbkey=Gene]
GTTTCGATCTTACCTGTGAAAGAACAAGACGCTGCCGACTGGCTAGCGCTGCGGAATCTTCTTTGGCTCG

Thank You

Answer: It's pretty straightfoward. What you'd probably find yourself caught out by is that biopython discards header info after the first space in the `sequence.id` field, so instead, you simply need to test the `sequence.description` field which retains it all.

i.e., if the first sequence is called `sequence` and is:

>lcl|KP771862.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]
ATGTTCAAACTTTTGAGTAAGTTATTGGTCTATTTGACCGCGTCTATCATGGCTATTGCG
AGTCCGCTCGCTTTTTCCGTAGATTCTAGCGGTGAGTATCCGACAGTCAGCGAAATTCCG
GTCGGGGAGGTCCGGCTTTA

Then `sequence.id` becomes: `lcl|KP771862.1_gene_1`, and `sequence.description` is `lcl|KP771862.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]`

Is just a matter of checking if "`blaVIM`" is `in` `sequence.description`.

>>> from Bio import SeqIO
>>> for rec in SeqIO.parse('bla.fasta', 'fasta'):
... if "blaVIM" in rec.description:
... print(rec.format('fasta'))
...
>lcl|KP771862.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]
ATGTTCAAACTTTTGAGTAAGTTATTGGTCTATTTGACCGCGTCTATCATGGCTATTGCG
AGTCCGCTCGCTTTTTCCGTAGATTCTAGCGGTGAGTATCCGACAGTCAGCGAAATTCCG
GTCGGGGAGGTCCGGCTTTA

>lcl|KP749829.1_gene_1 [gene=blaVIM] [location=1..801] [gbkey=Gene]
ATGTTCAAATTTTTGAGTAAGTTATTGGTCTATTTGACCGCGTCTATCATGGCTATTGCG
AGTCCGCTCGCTTT
Question: Hi I am trying to put together a plink ped file but I couldn't find how to depict missing maternal and paternal ids. You can give an argument if you have no parent ids in your file but I only miss some. For sex and phenotype you have missing data codes but what about the parental ids?
To put it differently how exactly do you depict the family relationships in the ped file format for plink?
Thanks

Answer: did you try to use '0' for the parental IDs ?
Question: I want to identify a list of accessible genes in a DNase-seq experiment. I downloaded the alignment file `ENCFF001CTZ.bam` from [Encode][1]. I performed peak calling by using macs1.4 and obtained a peaks.bed file containing coordinates of peaks. Next, I performed peak annotation by Homer using `annotatePeaks.pl`. The Final annotation file I obtained contains following columns:

"Peak_ID" "Chr" "Start" "End" "Strand" "Peak.Score" "Focus.Ratio.Region.Size" "Annotation" "Detailed.Annotation" "Distance.to.TSS" "Nearest.PromoterID" "Entrez.ID" "Nearest.Unigene" "Nearest.Refseq" "Nearest.Ensembl" "Gene.Name" "Gene.Alias" "Gene.Description" "Gene.Type"

A screenshot of Annotation_File from Homer is [here][2].

**QUESTION:** How I can finally obtain a list of genes which I can say are accessible? Do I need to only take genes for whose the column '*Annotation*' tells me that region is near to '*promoter-TSS*' or I can also take those which are '*Intergenic*' or so?

What information this (that a gene's intergenic region or intron is accessible) can provide?

[1]: https://www.encodeproject.org/experiments/ENCSR000EMT/
[2]: https://www.dropbox.com/s/0v1iflhcmbtu01l/Screenshot%20from%202015-12-08%2017%3A22%3A44.png?dl=0

Answer: Strictly speaking, any DNAseI HS site that overlaps a gene would be an "accessible gene".

But I'm guessing that you're looking for a transcribed gene(?); a gene that is accessible to PolII? If that's the case, you want the promoter-TSS overlaps.

DNaseI sites that sit elsewhere in the genome could be enhancers, but you would need more information to define these. Or they could be virtually anything else. I'm not aware of any studies that have tried to characterise the protein profiles associated with HS sites at non-promoter sites.
Question: I try to split fastq files into 'N' chunks using a simple CC program and htslib-C . It works fine:

./split2file -o TMP S1.R1.fq.gz S1.R2.fq.gz -n 10

but when I use a thread pool ( As far as I understand I can use one and only one thread pool for reading and writing the files ?)

./split2file -o TMP S1.R1.fq.gz S1.R2.fq.gz -n 10 -t 2
[LOG] closing "TMP.00001.R1.fastq.gz" N=100.
Segmentation fault

When I use threads I've got a SIGSEGV signal.

here is the gdb stacktrace:

```
#0 0x00007ffff718209b in pthread_cond_signal@@GLIBC_2.3.2 () from /lib64/libpthread.so.0
#1 0x00007ffff7f56770 in hts_tpool_dispatch3 () from /CONDAS/users/ifb-tools/miniconda3/envs/samtools-1.15.1/lib/libhts.so.3
#2 0x00007ffff7f11d42 in mt_queue () from /CONDAS/users/ifb-tools/miniconda3/envs/samtools-1.15.1/lib/libhts.so.3
#3 0x00007ffff7f148e5 in bgzf_flush () from /CONDAS/users/ifb-tools/miniconda3/envs/samtools-1.15.1/lib/libhts.so.3
#4 0x00007ffff7f14b76 in bgzf_close () from /CONDAS/users/ifb-tools/miniconda3/envs/samtools-1.15.1/lib/libhts.so.3
#5 0x0000555555557a3d in PairedEnd::~PairedEnd() ()
#6 0x0000555555556bc3 in main ()
```

here is the code:

(removed)

What Am i doing wrong ?

Answer: answering my question:

i was calling `hts_tpool_destroy(pool);` before closing the bgz output writers.
Question: I would like to use bedtools merge to collapse together all the features sharing a same `gene_id` in my bed file (which contains the annotation of various genes - the name column (4th one) also corresponds to the `gene_id`). Due to splicing, elements may be quite distant....
​
```
chr14 49894259 49895806 ENSMUST00000053290 0.000000 - mm10_ensGene exon . gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14 49894873 49894876 ENSMUST00000053290 0.000000 - mm10_ensGene stop_codon . gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14 49894876 49895800 ENSMUST00000053290 0.000000 - mm10_ensGene CDS 0 gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14 49895797 49895800 ENSMUST00000053290 0.000000 - mm10_ensGene start_codon . gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
chr14 49901908 49901941 ENSMUST00000053290 0.000000 - mm10_ensGene exon . gene_id "ENSMUST00000053290"; transcript_id "ENSMUST00000053290";
```

Answer: Hi- You don't say what to do with the collapsed features. This piece of code collapses the gene ids and counts how many features have been grouped together along with the min and max coordinates. Note it uses bedtools groupby not merge.

Sample data:

```
cat genes.bed
chr14 49894259 49895806 ENSMUST00000053290 0.000000 ...
chr14 49894873 49894876 ENSMUST00000053290 0.000000 ...
chr14 49894876 49895800 ENSMUST00000053291 0.000000 ...
chr14 49895797 49895800 ENSMUST00000053291 0.000000 ...
chr14 49901908 49901941 ENSMUST00000053291 0.000000 ...
```

Example output:

```
sort -k4,4 genes.bed \
| groupBy -g 1,4 -c 4,2,3 -o count,min,max \
| awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'

chr14 49894259 49895806 ENSMUST00000053290 2
chr14 49894876 49901941 ENSMUST00000053291 3
```
Question: I received genome-wide association (GWAS) data from a colleague who's supposedly done all the imputation and quality control according to the consortium's standards. Genotyping was Illumina 660, imputed to HapMap (3.2 million SNPs total).

The data came to me as a matrix of 11,000 samples (rows) and 3.2 million SNPs (columns). There's a header row for each SNP, and genotypes are coded as the number of minor alleles (or allele dosage for imputed SNPs).

Here's a few rows and columns to show you what it looks like:

rs1793851 rs9929479 rs929483
2.0 0 1
1.6 0 1
2.0 NA 0
2.0 0 1
1.6 0 0
2.0 1 NA
1.0 0 0
1.9 0 2

I've always used PLINK for GWAS data management, QC, and analysis because of its efficient data handling capabilities for GWAS data. However, this kind of data can't be imported directly into PLINK or converted into a pedigree format file. (PLINK does handle imputed data, and so does SNPTEST, but both of these require genotype probabilities and I only have the expected allele dosage).

I did write some R code to read in the data in chunks and run some simple summary and association statistics, but this is clunky and suboptimal for many reasons:

The dataset first has to be split up (I used a perl wrapper around UNIX/cut to do this). After splitting the dataset into several hundred files with all my samples and a subset of SNPs, computing sample-level measures (sample call rate, relatedness, ethnic outliers) is going to be a real coding nightmare.

Subsetting analyses is going to be difficult (not as easy as PLINK's --exclude, --include, --keep, --remove, --cluster, etc).

PLINK integrates SNP annotation info (in the map file) to your results. Joining QC and analysis results to genomic position, minor allele, etc, will require lots of SQL joins.

Ideally I don't want to rewrite software for GWAS data management, QC, and analysis. I've considered (1) analyzing only genotyped SNPS, or (2) rounding the allele dosage to the nearest integer so I can use PLINK, but both of these methods discard useful data.

Does anyone have any suggestions on how I should start to QC and analyze this data without re-inventing the wheel or rewriting PLINK? Any other software suggestions that could take this kind of data? Keep in mind, my dataset is nearly 100GB.

Thanks in advance.

Answer: Actually PLINK does now support dosage formats. a) You need to set format=1 and b) the dosage column names needs to "FID IID" format.

plink --dosage test.trdose format=1 --fam test.fam

Of course you need to transpose the mldose file (and rename "Al1" to "A1" and "Al2" to "A2"), as you would with the mlprob file. See http://pngu.mgh.harvard.edu/~purcell/plink/dosage.shtml for other option arguments.
Question: Can someone please help me count the number of reads in fastq file with 3' terminal bases? For example, I need to get the account of reads ending with A, T, G and C in any input.fastq

Answer: Heng Li's `seqtk` collaborating with awk (or here mawk which is faster):

seqtk seq -A in.fastq | mawk '$1 ~ ">" {print ">foo";next} {n=split($1,A,""); print A[n]}' | seqtk comp | mawk 'OFS="\t" {sumA+=$3; sumC+=$4; sumG+=$5; sumT+=$6} END {print "A:"sumA,"C:"sumC,"G:"sumG,"T:"sumT}'

It basically creates a new fasta file, all with the dummy name `>foo` (just that `seqtk` does not complain about missing headers), and the last base of the sequence, which is then sent to `seqtk comp`, summarizing the number of occurrences of ATCG in that file. Last part is simply adding each ATCG and printing it.

Still, I encourage you to try and find these kinds of solutions yourself. In my experiences, putting together these little scripts really help to make yourself better with the command line. The solutions are out there in the communities, you only have to pipe the individual pieces together to get your specific output.
Question: I have PLINK dosage files in the form pgen/psam/pvar. I would like to know how to do either/both of the following:

1. Convert the file set to bed/bim/fam files (hard-calls)

2. Read the allele dosage from pgen into R as a continuous variable

Thanks.

Answer: Convert to hard-calls: "plink2 --pfile ... --make-bed --out ..."

Convert to an R-readable form: "plink2 --pfile ... --export A-transpose --out ..."
Question: Hello everyone,
I was looking at some PCA plots taken from an experiment. These data don't cluster in their first two principal components and if I dig a little bit deeper int he analysis and look at other combinations I can see a clear clustering when comparing 5 (3% of variance) vs 6 (2% of variance).

What does that mean? We would have expected a big separation in the first 2 main principal components. Does it mean that the expected separation is minimal and is explained by the (3+2)%=5% of variance of the whole experiment?

Thanks in advance.

Answer: Could it be that there are technical confounders leading to more variability than the biological differences?
Question: Hi,

I'm trying to annotate genomic intervals with their corresponding genes in R.
My code:
```
db = TxDb.Hsapiens.UCSC.hg38.knownGene
bed = import.bed(bed_path)
genes = transcriptsByOverlaps(ranges = bed, x = db, maxgap=0, columns='gene_id')
```
I get for each interval the overlapping part to each gene.
I want to get the original interval and the corresponding genes.

For example if my BED file look like this:
```
chr1	152078793	156262976
```
I want to get:
```
chr1	152078793	156262976 gene1
chr1	152078793	156262976 gene2
chr1	152078793	156262976 gene3
```
and not:
```
chr1 152079557 152080903 gene1
chr1 152084141 152089064 gene2
```

Like bedtools intersect with the `-wo` parameter.

Thank you!

Answer: My bad. I misinterpreted. U can try below thing.

mergeByOverlaps(bed, genes)
Question: I have a combined fasta file with all my sequences. I want to print the ID lines and their DNA sequence that end in `.4`

So far I have `awk /'>*.4'/ {getline;print} combined.fasta` Which prints the sequences that I want, how do I get the ID lines as well?

Answer: grep '^>.*4$' -A 1 --no-group-separator in.fa
Question: What possible features can be extracted from nucleotide sequences for machine learning?
For example gc content, dinucleotide frequency etc.

Answer: There are multiple ways to compile your features:

1) ***Knowledge-based approach***: here you would only use a set of limited features that have a direct influence on your prediction/classification/learning task. Feature set will be limited, and you won't be able to add new knowledge to the field.
[See][1] an example where we used a subset of features that we assumed to have role in 3D domain swapping

2) ***Data-driven approach***: you can compile all available features that you can gather from your nucleotides (DNA or RNA?) and test them using rigorous feature selection method
[See][2] an example where we used the entire set of features with hybrid features (combining multiple features) to predict 3D domain swapping

3) ***Feature engineering/representation learning***: you can either of the above set and use [deep neural encoding][3] methods, here the algorithm would engineer the features ([NN, RBM, PCA, LSTM, etc.][4]). This approach is more applicable when you have large dataset(s) and not primarily looking for features contributing to your predictive model such as feature selection or biological inference.

- Count(s) of individual bases (ATGC - mono, bi, tri...)
- k-mer count (See previous answers)
- [physicochemical properties][5] of your sequences
- evolutionary scores (example [here][6])
- mutation/substitution scores (GERP, PhyloP, etc.)
- [Annotation][7]-based features (part of gene-structure (exon-intron), coding or [non-coding][8] etc.)

PS. Like one of the answers, it all depends on your prediction problem

[1]: http://www.ncbi.nlm.nih.gov/pubmed/20634983
[2]: http://www.ncbi.nlm.nih.gov/pubmed/21592079
[3]: https://t.co/YTs2xDLcai
[4]: https://en.wikipedia.org/wiki/Feature_learning
[5]: http://www.annualreviews.org/doi/abs/10.1146/annurev.bi.36.070167.002203
[6]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3145250/
[7]: http://www.ncbi.nlm.nih.gov/pubmed/26494363
[8]: http://www.nature.com/nmeth/journal/v11/n3/full/nmeth.2832.html
Question: I want to parse XML response obtained from bioproject DB using efetch module in Biopython.

Here is my code:

```py
from Bio import Entrez
Entrez.email = "myemail@company.org"
handle = Entrez.efetch(db="bioproject", id="55465", rettype='gb',retmode="xml")
records = Entrez.parse(handle)
for record in records:
print record
```

but this gives the following error:

```py
Bio.Entrez.Parser.ValidationError: Failed to find tag 'RecordSet' in the DTD. To skip all tags that are not represented in the DTD, please call Bio.Entrez.read or Bio.Entrez.parse with validate=False.
```

instead if I try this, it works but it gives the XML lines as it is (no parsing)

```py
handle = Entrez.efetch(db="bioproject", id="55465", rettype='gb',retmode="xml")
readlines = handle.readlines()
for line in readlines:
print line
```

Can anyone please guide as what is the right way to parse the XML response in this case?

Answer: This seems to be a problem with the Bioproject database not having a DTD: http://lists.open-bio.org/pipermail/biopython/2015-May/015632.html

That thread has a work-around, but you might also want to try the development version of Biopython (https://github.com/biopython/biopython) to see if the new work Peter mentions will "just work".

**EDIT**

I'm not sure it will make any difference, but 'gb' is not a rettype available for the Bioproject database:

http://www.ncbi.nlm.nih.gov/books/NBK25499/table/chapter4.T._valid_values_of__retmode_and/?report=objectonly
Question: I need clinical data of Glioblastoma patients from TCGA including Age at initial pathologic diagnosis, Date of death, Last follow up date, Ethnicity, Gender, Histological type, History of neoadjuvant treatment, Method of initial pathologic diagnosis, Karnofsky performance score, Performance status scale timing, Person neoplasm cancer status, Prior glioma, Race. How to download these data from TCGA portal giving patient or sample ID's? If anyone have this dataset, plz share as soon as possible.

Answer: You can do that in R using the TCGAbiolinks package (with which you can also download the actual gene expression data as [follows][1]) and follow the instructions from the [Bioconductor webpage][2]. My script would be as follows:

library("TCGAbiolinks")

DIRPREFIX <- "/PATH/"

## Download raw counts for indicated sets
## Example of TCGA.data.sets: c("TCGA-GBM")
for (id in TCGA.data.sets){
clin.query <- GDCquery(project = id, data.category = "Clinical")
json <- tryCatch(GDCdownload(clin.query),
error = function(e) {GDCdownload(clin.query, method = "client")})
## the following types of clinical data can be obtained (just pick whichever you need)
clinical.drug <- GDCprepare_clinic(clin.query, clinical.info = "drug")
clinical.admin <- GDCprepare_clinic(clin.query, clinical.info = "admin")
clinical.radiation <- GDCprepare_clinic(clin.query, clinical.info = "radiation")
clinical.patient <- GDCprepare_clinic(clin.query, clinical.info = "patient")
clinical.stage_event <- GDCprepare_clinic(clin.query, clinical.info = "stage_event")
clinical.new.tumor.event <- GDCprepare_clinic(clin.query, clinical.info = "new_tumor_event")
clinical.followup <- GDCprepare_clinic(clin.query, clinical.info = "follow_up")
clinical.index <- GDCquery_clinic(id, type = "clinical")
clinical.biospecimen <- GDCquery_clinic(id, type = "biospecimen")
save(WHAT_YOU_NEED, file = file.path(DIRPREFIX, paste0(id, "_clinical.Rdata")), compress = "xz")
}

[1]: https://www.biostars.org/p/211956/#215791
[2]: http://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html