cannin/biostars_qa huggingface dataset:
Question: Good day all,

So I have been trying to retrieve RNA-seq, miRNA-seq and protein expression data from TCGA' GDC data portal. Unfortunately, I have been having some issues. I have read a number of previous posts here in biostars and other forums, but still having issues using the new GDC portal

I have figured out how to get the Primary Tumor and Solid normal tissue RNA-seq and miRNA-seq datasets I need i.e. TCGA barcode ids with 01 or 11 code respectively - as defined by
https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes

For example TCGA-COAD project has 456 RNA-seq files and 444 miRNA-seq files

However, I am unable to find the corresponding protein expression datasets - especially Solid normal tissue datasets(11 code or any other). I have looked at TCPA - https://tcpaportal.org/tcpa/
and only found Primary tumor protein expression (01 code) for the TCGA projects I need. I have also looked through the legacy GDC portal but was unable to find any way to match the files names to the TCGA barcodes.
In the new GDC portal I can put files in my cart and download the "sample sheet" which gives me the mapping of actual filenames, MD5 hashcodes, TCGA barcodes, etc. I honestly was under the impression that GDC would have the protein expression data as listed in this link - https://cancergenome.nih.gov/abouttcga/aboutdata/datalevelstypes

Am I missing something here? are my assumptions wrong on that Solid normal tissue datasets(11 code) are included? Amy help would be greatly appreciated.

Answer: **Edit April 5th, 2020:**

Since posting this answer, the NIH has set up the Proteomic Data Commons

-------------
----

As far as I am aware, protein expression profiling was not a primary goal of TCGA. So, some cohorts may not even have this type of data. Thus, if you cannot find any normals in your data of interest, that may be because no normals were ever profiled.

I would also take a look a cBioPortal, which has RPPA (reverse phase protein array) data for different datasets, including colon, I see.

Kevin
Question: I have a bgzip vcf file and build tabix index already, but I am wondering if tabix can query VCF file by ID, just like:

```
tabix my.vcf.gz id='rs527352437'
```

or is there another way to query dbSNP by ID quickly?

Answer: Here's an approach that uses `hg38` SNP data from UCSC, an `rs*` ID of interest, and `pts-line-bisect` to do a binary search.

Your VCF-formatted data will be in different columns, but the principle here is the same — you permute the ID into the first column of a text file. You sort that text file on the first column. You do a binary search on the sorted text file.

Binary searches on sorted data are fast — *log n* fast — resulting in searches that take fractions of a second where they would otherwise take minutes or longer.

1) Get SNPs and write them into a text file sorted by SNP ID. For `hg38`, for instance:

$ LC_ALL=C
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/snp147.txt.gz \
| gunzip -c \
| awk -v OFS="\t" '{ print $5,$2,$3,($3+1) }' \
| sort -k1,1 \
> hg38.snp147.sortedByName.txt

2) Install [pts-line-bisect][1]:

$ git clone https://github.com/pts/pts-line-bisect.git
$ cd pts-line-bisect && make && cd ..

3) Run a query:

$ rs_of_interest=rs2814778
$ ./pts-line-bisect/pts_lbsearch -p hg38.snp147.sortedByName.txt ${rs_of_interest} \
| head -1 \
| cut -f2- \
> result.bed

In your case, you could modify your VCF file to permute the columns, so that the ID is in the first column of all columns. You permute and then sort as in step 1. This is the file you query with `pts-line-bisect`. Whatever result you get — if any — you permute back into VCF field ordering, likely using `awk` or similar.

With some modifications, it should be possible to do a binary search on a sorted, compressed text file. An important consideration is that the sort order is stable, so that compressed blocks provide a consistent or "deterministic" direction to go in finding the right compressed block for the data being queried.

[1]: https://github.com/pts/pts-line-bisect
Question: I'm trying to use WGCNA for weighted network analysis. I'm getting trouble how to decide this parameter. It's clear how to define and why this has been defined.

#Choosing a soft-threshold to fit a scale-free topology to the network

powers = c(c(1:10), seq(from = 12, to=20, by=2));

What is the meaning of this powers and how to decide this value?

Answer: In WGCNA, power parameter is needed to reduce the spurious correlations in the data. To select a power from `pickSoftThreshold` function, choose reasonably high `R^2` (column 2), usually higher than say .85 and negative slope around -1 (column 4) to get an approximately scale-free network. Power is helpful to easily differentiate strong and weak correlations.
Question: Below is a simple example of gff3 file:

1	T1	gene	3631	4605	.	+	.	ID=ATNG01010
1	T1	mRNA	3631	4605	.	+	.	ID=ATNG01010.1;Parent=ATNG01010
1	T1	exon	3631	3913	.	+	.	ID=ATNG01010:exon:1;Parent=ATNG01010.1
1	T1	CDS	3860	3913	.	+	0	ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1	T1	exon	3996	4276	.	+	.	ID=ATNG01010:exon:2;Parent=ATNG01010.1
1	T1	CDS	3996	4260	.	+	2	ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1	T1	exon	4486	4605	.	+	.	ID=ATNG01010:exon:3;Parent=ATNG01010.1

My question is: if we found another coding sequence (encode a different protein) range from 3752 to 3904, how should the gff3 file look like? It seems to me that the gff3 file can only allow one protein-coding gene per mRNA. If not, could anyone show me one example? Thank you!

Answer: In principle, you can define an arbitrary number of CDS per mRNA. The Parent attribute of each CDS indicates to which mRNA it belongs. If your CDS feature spans multiple lines (discontinuous features) it must have an ID to indicate lines that collectively represent the CDS. In fact, your example already contains two different protein coding sequences for mRNA 'ATNG01010.1', namely 'ATNG01010:CDS:1' and 'ATNG01010:CDS:2'. You could add a third one using the same pattern.

GFF Specification: https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md
Question: I'm lost at the NCBI FTP server. Which file contains all chromosomes of GRCh37.p13 in one fasta? The program that I want to use takes one fasta file as the reference sequence. At ftp://ftp.ncbi.nih.gov/genomes/Homo_sapiens/ARCHIVE/ I don't see GRCh37.p13 and I don't see a README file .. I'm a little bit stuck here.

Answer: The NCBI record for GRCh37.p13 (GCF_000001405.25) can be found here: https://www.ncbi.nlm.nih.gov/search/?expt_dict=b&term=grch37.p13.

Click on *Download* and then select *Genomic FASTA* and *RefSeq*:

The tar file contains the complete NCBI version of the genome release in FASTA format.

Kevin

**NB - FTP site is here: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Homo_sapiens/all_assembly_versions/**
Question: I can't figure out how to use the "-match" syntax even after reading all the documentation I could find. I get these errors:

```
$ cat xml.txt | xtract -pattern Gene-commentary -match Gene-commentary_type:1

Unrecognized argum`ent '-match'
No -element before 'Gene`-commentary_type:1'

$ cat xml.txt | xtract **-element** Gene-commentary -match Gene-commentary_type:1

Unrecognized argument '-match'
No -element before 'Gene-commentary_type:1'
```

**What am I doing wrong?**

---

What I am trying to do is pull the accession of the reference sequence and the coordinates for the region for a given entry in NCBI Gene (see https://www.biostars.org/p/122522/) so that I can run `efetch -format FASTA -seqstart -seqend` and get the appropriate results.

I could parse the XML in python to do it, but it really seems like I should be able to do this in "one line" using entrez direct if only I could get `-match` to work :/

Here is what the XML looks like:

Say I have a gene record in XML

epost -db gene -id 672 | efetch -format xml > xml.txt

According to the outline,

```
cat xml.txt | xtract -outline

1
Reference assembly
RefSeqGene
NG_005905
2

92500
173688
```

I have read:

https://www.biostars.org/p/94539/

http://www.ncbi.nlm.nih.gov/books/NBK179288/ (I followed these instructions to install it, so `which epost` returns `~/edirect/epost`)

http://www.ncbi.nlm.nih.gov/news/02-06-2014-entrez-direct-released/?campaign=facebook-02072014

http://elane.stanford.edu/laneconnex/public/media/documents/EntrezDirect.pdf

Answer: Using a good old xslt stylesheet:

```

```

run:

xsltproc --novalid transform.xsl "http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=gene&id=672&retmode=xml"

output:

```
(Reference GRCh38 Primary Assembly)NC_000017:43044294-43125482
(Reference assembly)NG_005905:92500-173688
(Alternate CHM1_1.1)NC_018928:41431850-41513017
(Alternate HuRef)AC_000149:36962662-37043808
```
Question: Hi all,

Plink has the function "Linkage disequilibrium based SNP pruning" which is `--indep 50 5 2` where the 2 stands for the vif threshold (VIF is 1/(1-R^2)) which means in this case r2 = 0,50.

*So Linkage disequilibrium* is the non-random *association* of alleles. I'm a bit strugling what for example an r2 of 0.5 means and how plink calculates this. Does 0.5 just mean a correlation of 0.5 between two snps? Can anyone explain to me a bit more what this 0.5 actually mean in real numbers. For example when I have 100 samples how many snps need to be in perfect LD to reach a r2 of 0.5?

Answer: It's the squared correlation coefficient between the 0/1/2 allele counts. I.e. ((Cov(marker 1 allele counts, marker 2 allele counts))^2) / (Var(marker 1 allele counts) * Var(marker 2 allele counts)).

The "number of SNPs in perfect LD" required depends on the exact distribution of allele counts for each marker; for example, if both markers have empirical MAF 0.01, even 99 samples in "perfect LD" is not enough to guarantee r^2 >= 0.5.
Question: I have RNA-seq samples that are highly contaminated with Illumina PCR primer 2.01, which program is best to use to remove those?

Answer: I'd recommend bbduck from the BBMap suite. It has included a lot of adapters and primer sequences, is very fast, and you can trim next to adapters also low-quality tails. Additionally, it is keeping your paired-end data in the correct order (if one read is removed, the other is removed too).
Question: Hi,

I have a huge list of p-values ranging from E-05 to E-324. I want to use a toll (Cropper) but it requires decimal numbers as input.
How can I transform scientific notations to decimals, all at once?

Answer: I'd first question if your program really doesn't accept scientific notation, if so that creates problems for working with these values and you should consider using a properly implemented program. The following example shows what I mean: In principle, the following bash code does the conversion just fine:

while read -r line;
do
printf "%.325f\n" $line ; ## change format to %.308f as minimum sensible representation
done < testfile.txt

You need to decide for some sensible cutoff of precision, I am using 325 just to show the problem:

Input (your range):

1E-5
1E-100
1E-324

Output:

0.00001000000000000000081803053914031309545862313825637102127075195312500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000000000000000199918998026028836196477607885341594201826030059365956992555434676176762886132929895827460748109118507985282705397496540222684360419612636083562831412787179427249289424690806658916305930004345786023014502507945

0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Note: 1E-324 is effectively zero in this notation but maybe smaller than the smallest positive non-zero double precision value. On 64bit platforms that is 2e-308 **

Of course you could do some trickery by having variable precision floating point numbers, but you rather don't want that.

**
checked with the c program:

#include
#include
int main()
{
printf("%.e\n", DBL_MIN);
}
Question: Hello,

I am working with the Biomark HD Fluidigm high throughput qPCR transcriptomics datasets (consists of 240 target genes and 8 housekeeping genes). The dataset consists of 250 samples processed in 5 different batches (50 samples/batch). Dataset of each batch was normalized to calculate `Delta Ct` (Difference between the Target gene and Geometric mean of 8 housekeeping gene (i.e GM_HKgs)). Then, transformed to `- Delta Ct` (negative delta Ct) because we want to apply downstream analyses like ANOVA to log-scale normalized relative expression values, as just using `-1 * Delta Ct` as readout is quite valid too. I understand from various posts `removebatcheffects` (`limma` package) would be best to use. My question is what data type should be used as input, Delta Ct or Negative Delta Ct (-1*Delta Ct).

Answer: You could use negative delta ct, which would be fine, or see https://support.bioconductor.org/p/67032/ for a more sophisticated normalization method that we use for qPCR.

You don't necessarily need `removeBatchEffect` however. The limma DE/anova analysis of qPCR ct values does not use `removeBatchEffects`.
Question: I am working with a program that generates artificial FASTQ files from a given reference genome and it allows for read length customization. I am trying to produce several FASTQ files with different read lengths. Hence, I want to know what read lengths are produced by real sequencers such as those from Illumina. I found [this][1] site which states the maximum read length of some sequencers.

However, when stating, for example, that the maximum read length is 2 x 150 bp, is that supposed to mean that each end of a contig will be 150 bp long? Also, if that is the case, does this mean that the sequencer could be set up to produce read lengths equal to any number of base pairs less than 150? i.e. could the sequencer produce read lengths of 2 x 1 bp, 2 x 2 bp, 2 x 3 bp,...., 2 x 148 bp, 2 x 149 bp? Lastly, which one of those sequencers is the most commonly used for genomes which are about as long as, say, that of E. Coli?

[1]: http://www.illumina.com/systems/sequencing.ilmn

Answer: The read lengths (up to maximum) are somewhat configurable. However, there are some very commonly used lengths that have predominated at different times. In earlier days you would see a lot of 2x36, 2x50, and 2x76 from the Illumina GA or GAII. With the HiSeq 2000/2500 it was most common to see 2x100 or 2x150 (you also often see 2x101 or 2x151). The MiSeqs often produce 2x250 or 2x300. The HiSeq 4000, HiSeq X, and NovaSeq 5000/6000 seem to most often be run at 2x150. I would probably go with 2x100 for your simulations.
Question: Hi,

Could provide the faster way to filter this data :

chr1	43	1000	gene_name=boby	gene_type=trucA
chr2	44	1000	gene_name=natt	gene_type=trucB
chr3	45	1000	gene_name=alurika	gene_type=trucC

To :

chr1 43 1000 boby trucA
chr1 44 1000 natt trucB
chr1 45 1000 alurika trucC

CORRECTION :
Original text data looks like this :

chr1 43 1000 TEST gene_name=boby;gene_type=trucA;foo=34
chr2 44 1000 TRUC gene_name=natt;gene_type=trucB;foo=34
chr3 45 1000 PASS gene_name=alurika;gene_type=trucC;foo=34

Answer: Thanks for your reply !

But I just discover right now that awk support regexp for the Fieldseperator.
So this works too :

cat test | awk 'BEGIN{FS="\t|="} {print $1,$2,$3,$5,$7}'
Question: Hi,

Will always R2 file contain sense transcript and R1 reverse complement of transcript?

If R2 file contains sense transcript, I will have to run HTSeq.scripts.count with "reverse" setting to quantify sense transcription. Quite confusing, isn't it?

Library creation kit -> E7420S NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina®

Thanks, Bernardo

P.S. Here is the pipeline:

set -ue

#################################################

#BWA test aligner with paired-end data

#################################################

# Get the genome file from the command line

genome_file=$1

# Get the fastq file from the command line

fastq_file_R1=$2

# Get the fastq file from the command line

fastq_file_R2=$3

# Get the fastq file from the command line

fastq_file_R3=$4

# Get the fastq file from the command line

fastq_file_R4=$5

#get gff

GFF=$6

#BWA default settings

bwa index $genome_file

#BWA input

bwa mem -t 8 $genome_file $fastq_file_R1 $fastq_file_R2 | gzip -3 > P_S1_L001_aln-pe.sam.gz

bwa mem -t 8 $genome_file $fastq_file_R3 $fastq_file_R4 | gzip -3 > V_S1_L001_aln-pe.sam.gz

################################################################

#Flagstat

################################################################

#Convert .sam to .bam to input to Flagstat

samtools view -b -S -o P_S1_L001_aln-pe.bam P_S1_L001_aln-pe.sam.gz

samtools flagstat P_S1_L001_aln-pe.bam

samtools view -b -S -o V_S1_L001_aln-pe.bam V_S1_L001_aln-pe.sam.gz

samtools flagstat V_S1_L001_aln-pe.bam

################################################################

#Count reads mapped with htseq-count

################################################################

samtools sort -n V_S1_L001_aln-pe.bam invivo.sorted

python -m HTSeq.scripts.count -m intersection-nonempty -f bam -a 10 -t mRNA -i Parent -r name -s reverse invivo.sorted.bam $GFF | awk 'n>=5 { print a[n%5] } { a[n++%5]=$0 }' > invivo_R.counts

samtools sort -n P_S1_L001_aln-pe.bam plate.sorted

python -m HTSeq.scripts.count -m intersection-nonempty -f bam -a 10 -t mRNA -i Parent -r name -s reverse plate.sorted.bam $GFF | awk 'n>=5 { print a[n%5] } { a[n++%5]=$0 }' > plate_R.counts

################################################################

#IGV

################################################################

#samtools sort -n sorts by name, not coordinate.. index requires sort by coordinate..

samtools sort V_S1_L001_aln-pe.bam V_S1_L001_aln-pe.bam.sorted

samtools index V_S1_L001_aln-pe.bam.sorted.bam

samtools sort P_S1_L001_aln-pe.bam P_S1_L001_aln-pe.bam.sorted

samtools index P_S1_L001_aln-pe.bam.sorted.bam

Answer: I believe the "how" is answered on this website (1st Figure under product information tab):

https://www.neb.com/products/e7420-nebnext-ultra-directional-rna-library-prep-kit-for-illumina

Seems like it is a dUTP protocol, which would mean you should be using "--library-type fr-firststrand" in TopHat:

http://ccb.jhu.edu/software/tophat/manual.shtml#toph

http://ccb.jhu.edu/software/tophat/faq.shtml#library_type
Question: I am using Pysam to process a bam file containing paired reads. However, I get the error

AttributeError: 'csamtools.AlignedRead' object has no attribute 'rnext'

when I try to extract the chromosome that the mate of a read is aligned to. What am I doing wrong?

for row in samfile:
print(samfile.getrname(row.rname)+' '+samfile.getrname(row.rnext))

It gives the read's reference properly, however it gives an error for the mate's reference.

Answer: Hi,

Looks like your pysam installation is too old and doesn't support rnext (see also QA3.2 [here](https://github.com/BSSeeker/BSseeker2)). Which version are you using? According to the logs rnext was introduced end of 2011, i.e. any version >= 0.6 should work.

Andreas
Question: I am a clinician, and have no particular expertise in genomics. I'm confused by much of what I see on dbSNP. Any help would be greatly appreciated! Thank you in advance.

I need some help understanding HGVS Names!

I see SNP's referred to with conventions like **A66G** for MTRR ([rs1801394][1]) and **A1298C** for MTHFR ([rs1801131][2]).

When I look under **HGVS names** for rs1801394 I see: NM_002454.2:c.66A>G, which l assume equals "A66G", but there are many other "names" there. What are these names? Are they all equivalent?

When I look under **HGVS names** for rs1801131 I do not see anything resembling "A1298C". Does this mean I have the wrong rs# for this SNP?

According to [SNPedia][3], it's the correct rs#: "rs1801131 is a SNP in the MTHFR gene, representing an A>C mutation at mRNA position 1298".

[1]: http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1801394
[2]: http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1801131
[3]: http://www.snpedia.com/index.php/Rs1801131

Answer: I'll try to explain briefly, but it's a lot and it's wonderful - one of the rare parts where we have great standards.Check it out: http://www.hgvs.org/mutnomen/

The usual format is:

:.

`NM_002454.2:c.66A>G`is a cDNA change of Adenine to Guanine at position `66` in the cDNA ref sequence `NM_002454` version 2.

`A66G` most commonly refers an amino acid change of Alanine to Glycine in some protein with some identifier (should ideally be `NP_xxxxxxx:p.A66G`). In your case though, I think you might be referring to a `66A>G` nucleotide change as `A66G`, which is wrong.

Ref sequences can be `g.` (genomic), `c.` (cDNA), `r.` (RNA) or `p.` (protein)

Changes can be single amino acid variants (`p.Ala66Gly`) single nucleotide variants (`A>G`), deletions(`c.10delA`), insertions (`c.10_11insT`), duplications (`c.10dupA`), indels(`c.10_11delACinsTGA`) etc.

It is so well standardized that programs can parse and pick individual components for us to use.
Question: I have a file in GFF format and I need to convert it to BED format. What do I do?

Answer: Both formats are tab delimited text files used to represent DNA features in genomes. The order of columns between the two are different, there are also columns that correspond to attributes missing from one or the other format. Nonetheless the most important difference between the two is the coordinate systems that they assume.

The BED format developed at UCSC uses a zero based indexing and an open end interval whereas the GFF format developed at Sanger assumes a 1 based coordinate system that includes both start and end coordinates. Therefore

The [0,100] interval in BED format corresponds to [1,100] in GFF format and both are 100 base long. That the first element in BED format will be have the index of 0 where the last 100th element will have the index of 99! Whereas in GFF the first element will have the index of 1 and the last element will have the index of 100.

To convert between the two you may use Galaxy and select the section called Select Formats that will list various transformation options.
Question: Quick question: is there out there a .bed (or similar) files that span over all coding regions in hg38 coordinates? I need to analyze some WGS samples

Thank you very much in advance for any help!

Answer: wget -O - "http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.annotation.gtf.gz" |\
gunzip -c | grep 'transcript_type "protein_coding"' |\
awk '($3=="exon") {printf("%s\t%s\t%s\n",$1,int($4)-1,$5);}' |\
sort -T . -t $'\t' -k1,1 -k2,2n | bedtools merge

?
Question: The peaks of a particular Transcription factor CHIP-seq as identified by multiple peak calling algorithms e.g. MACS2 peakzilla e.t.c identify peaks of varying sizes or regions from 50 bp up to 1000 bp or more. In that case I was wondering that if the average size of transcription factor binding site is 10 bp based on available data in OReganno e.t.c (as indicated by [this post][1]) would it be safe to assume that for a particular peak the TFBS location is in its center or not. If not then what would be a good way to go about doing it and if such a method already exist or if someone has already done it.

[1]: https://www.biostars.org/p/64854/

Answer: No, it would not be safe to assume that the peak is at the center. If you know the motif for your protein, you can simply search for that motif in your ChIP-seq peaks. If you do not know the motif, you can use any number of software packages that will look for enriched sequences in the many peaks that you have to try to define what the motif is. In both cases, after you have the motif, then you "know" the location.