cannin/biostars_qa huggingface dataset:
Question: I tried to use BWA MEM to map reads from an interleaved FASTQ.

fastq="all.fastq"
fasta="/share/PI/apps/bcbio/genomes/Hsapiens/GRCh37/seq/GRCh37.fa"
bwa="/share/PI/apps/bcbio/anaconda/bin/bwa"
nThreads="12"

#Run BWA MEM
#IMPORTANT: NEED -p since "$fastq" is an interleaved fastq
readGroup="@RG\tID:CHM1\tSM:CHM1\tPL:Illumina"
sam="CHM1.sam"
"$bwa" mem -R "$readGroup" -t "$nThreads" -p "$fasta" "$fastq" -o "$sam"

(The FASTQs are CHM1; I used `prefetch` to fetch `.sra` files from three different runs from NCBI, then used `fastq-dump` to convert the SRAs to FASTQs, then `cat`ed them all together into one FASTQ.)

The SAM is 515Gb but has no obvious problems. `samtools quickcheck` says it's valid. But when I run GATK4's `FixMateInformation` or `ValidateSamFile`, I get output like this

ERROR: Record 1, Read name ######################################################################################################################################################################################################, Zero-length read without FZ, CS or CQ tag
WARNING: Record 1, Read name ######################################################################################################################################################################################################, QUAL field is set to * (unspecified quality scores), this is allowed by the SAM specification but many tools expect reads to include qualities
ERROR: Record 421522661, Read name ######################################################################################################################################################################################################, Zero-length read without FZ, CS or CQ tag

There may be even more errors, but this is what I got after two hours.

I can, in fact, see that the first line in the SAM file is

###################################################################################################################################################################################################### 4 * 0 0 * * 0 0 * * AS:i:0 XS:i:0 RG:Z:CHM1

Is this SAM really invalid? Or is there something I need to do so GATK4 will accept it?

Answer: `fastq-dump` will create file names and output reads to those files, so don't use redirection to stdout (`>`) - if you want to use it, then you have to use `-Z`:

-Z|--stdout Output to stdout, all split data become
joined into single stream

There is a working example of fastq-dump for the very same dataset you are using here:

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

fastq-dump --origfmt -I --split-files --gzip SRR1514950.sra

This will output two files one for R1 and other for R2, so you will need to modify your bwa mem command.
Question: Hi all. I'm working on Python script that find ORF for a sequence, within a function of code. I hope this is well explained enough to follow:

It works so far if there is a start and stop codon found, or just stop. However if there is a start but no stop codon I run into a problem with one of my sequences, the return from the function is "None".

Here's what I have thus far (might be indent issues from pasting):

def find_orf(sequence, gb):
start_pos = sequence.find('GCCGCCACCATG')
print "START " + str(start_pos)
if start_pos >= 0:
s_to_ATG = int(start_pos) + 9
start = sequence[s_to_ATG:]
for i in xrange(0, len(start), 3):
stops =["TAA", "TGA", "TAG"]
codon = start[i:i+3]
## start code block with issue
if codon in stops:
orf = start[:i+3]
else:
orf = start
## end code block with issue
return orf, start_pos
elif start_pos < 0:
stop_pos = sequence.find(str(gb[-12:]))
begin_to_stop = int(stop_pos) + 12
return sequence[:begin_to_stop], start_pos

else:
print "Error: There is no open-reading frame for this sequence!"

I've highlighted what is giving me issues in the code. It seems to always go to the "else" statement. The first sequence I have, I know for sure it has a START and STOP. The second sequence does not have a STOP, but has a START. So, I want it to print START to the end of the sequence, which in the code is the variable "start".

If I change this IF statement in the code to this:

for i in xrange(0, len(start), 3):
stops =["TAA", "TGA", "TAG"]
codon = start[i:i+3]
if codon in stops: return start[:i+3], start_pos

it prints the first sequence start to stop correctly, and the second is return as "Nonetype", because it has no STOP. It may be something simple that I'm overlooking, but I was hoping someone could see something wrong with the first code example, so that if there is a start and stop, it will print it correctly, and if there is no stop will print start to the end. (The second part of the code works where if there is no start, it prints from the beginning to the stop)

All help is appreciated.

Answer: You just need to add a break:

def find_orf(sequence, gb):
start_pos = sequence.find('GCCGCCACCATG')
print "START " + str(start_pos)
if start_pos >= 0:
s_to_ATG = int(start_pos) + 9
start = sequence[s_to_ATG:]
for i in xrange(0, len(start), 3):
stops =["TAA", "TGA", "TAG"]
codon = start[i:i+3]
if codon in stops:
orf = start[:i+3]
break ##Added a break statement here
else:
orf = start
return orf, start_pos
elif start_pos < 0:
stop_pos = sequence.find(str(gb[-12:]))
begin_to_stop = int(stop_pos) + 12
return sequence[:begin_to_stop], start_pos

else:
print "Error: There is no open-reading frame for this sequence!"

Without the break, you'll always iterate to the end of the sequence.
Question: I'm working on microarray data of GEO. Many data are in .txt format. how can I use them for analysis. Is there any way to convert .txt format to .CEL format?

Answer: To answer your question, no, there is not a way to convert .txt to .cel.
Question: Hi

I'm having a problem with a whole-exome-sequenced dataset consisting of about 400 human subjects, 200 cases with a certain disease, and 200 controls without. The dataset has been through a rigorous quality control (standardised QC in plink with HWE, IBD, missingness, sex-check ++, along with HapMap population stratification and Eigenstrat/PCA-analysis). I´m using plink to do a basic association-analysis for all variants between cases and controls, and while the resulting QQ-plot for the common (MAF > 0.01) variants is OK, the plot for the rare (MAF < 0.01) variants is less so. Below are the three QQ-plots for all, common and rare variants along with lambda-values:

- QQ-plot all variants, lambda 1.83
http://postimg.org/image/oyr53wchr/
- QQ-plot common variants, lambda 1.04
http://postimg.org/image/6lqjtc20v/
- QQ-plot rare variants, lambda 2.43
http://postimg.org/image/ic4haputb/

The main problem seems to be the positive deviation (observed > expected) of the rare variants in the first part of the plot, causing the lambda to be very big, both for the QQ-plot for rare and all variants. I am wondering what could be the cause of this behaviour for the rare variants, and also what the implications this has for the prospects of doing analysis on rare and common variants together.

I would be grateful if anybody has any experience in these matters and could provide some input.

Many thanks.

Answer: I would generate p-values using an exact test, and see if there is still inflation, personally.
Question: Hello,

I want to filter my non-variant positions for mapping quality. I used mpileup to output the mapping qualities across all sites. To do the filtering, can I take an arithmetic mean of the mapping qualities for each read that the base belong to?

For example:

NW_008793873.1 13 G 2 .^S. AB HS

Position 13 is covered by two reads with mapping qualities **H** and **S**. Would **(39 + 50)/2 = 44.5** be correct?

I ask this because in connection to mapping qualities, usually root mean square is mentioned so I was wondering what would the correct approach be in this case?

Thank you!

Answer: These are probabilities of mismapping on a PHRED scale. For the first one, the probability of mismapping is:

(10^(-(39/10)) = 0.0001258925

For the second it is:

10^(-(50/10)) = 1e-05

So on average, your probability of mismapping is:

(0.0001258925+1e-05)/2 = 6.794627e-05

On a PHRED scale it is:

-10*log10(6.794625e-05) = 41.67835
Question: Hello guys. I met some problems when I tried to follow the tutorial of Nanopolish in https://nanopolish.readthedocs.io/en/latest/quickstart_eventalign.html

One of the steps is to use samtools to sort alignments. Here I got a reads.bam file, and try to sort it into reads.sorted.bam file. My environment is Utubun and the version of samtools is 1.9 (the most recent one).

I run the following command

samtools sort -o reads.sorted.bam reads.bam

Error: samtools sort: couldn't allocate memory for bam_mem

Folowing the suggetion from https://github.com/samtools/samtools/issues/831
Then I try to run

samtools sort -m GALAXY_MEMORY_MB / GALAXY_SLOTS -o reads.sorted.bam reads.bam

Output information:

Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-n Sort by read name
-t TAG Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
-o FILE Write final output to FILE rather than standard output
-T PREFIX Write temporary files to PREFIX.nnnn.bam
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0]

For me, now it seems no error occur... However, I could not find the reads.sorted.bam in the current directory... so when I run the next step,

samtools index reads.sorted.bam

Expected error occurs: samtools index: failed to open "reads.sorted.bam"; No such file or directory.

I really dont understand what happened. Can anyone help me to find out why I could not create the sorted bam file? Thank you.

Answer: > Error: samtools sort: couldn't allocate memory for bam_mem

use option -m `(-m INT Set maximum memory per thread; suffix K/M/G recognized [768M])`

to reduce the memory. For example : '-m 10M' but it looks like your computer hasn't enough memory...
Question: hi,

I have a long list of different miRNAs like below

hsa-miR-7641
bta-miR-2904
hhi-miR-7641
hsa-miR-4454
bta-miR-2478
efu-let-7c
mmu-miR-6240
hsa-miR-7704

is there any way to avoid searching one by one in mirbase and copy pasting their mature sequences and taking their sequences at the same time?

Answer: Hi ,

You can use [faSomeRecords][1] to extract multiple fasta records

Usage:

faSomeRecords mature.fasta your-list-file result.fa

in the `grep -f/-F` be careful of space and blank lines. you can check them with [Notepad++][2] in windows.

Take care.

[1]: http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faSomeRecords
[2]: https://notepad-plus-plus.org/download/v7.html
Question: Hi all,

Sorry, if you find the question is basic. Could you please let me know why "A/T" and "C/G" called ambiguous SNP? why if such a SNP, say A/T, locate on one strand, we cannot T/A is on the opposite strand?

Thank you

Answer: So if you know that a SNP is A/T on the forward strand it is T/A on the reverse strand. That's fine. You can check the reference genome and it's clear which is the reference allele and which is the alternative.

The ambiguity is when someone says something like "A is the minor allele" or "A is associated with this phenotype" or "A is the ancestral allele" and they are not clear about which strand they are talking about. Then you don't know if they mean forward A and T is the minor/phenotype-associated/ancestral, or if forward T reverse A is the minor/phenotype-associated/ancestral.
Question: Hi all,

I am using Ensembl VEP (command line) to annotate a VCF I have. I am specifically looking for gnomAD allele frequencies, which is fairly straight forward to do, technically speaking. However, the data looks off in some cases.

For example, when I pass in:

10	69408929	COSM3751912	A	T	.	.	GENE=TACR2;STRAND=-;CDS=c.734T>A;AA=p.M245K

I get the VEP output:

COSM3751912	10:69408929	T	ENSG00000075073	ENST00000373306	Transcript	missense_variant	1278	734	245	M/K	aTg/aAg	rs55953810	MODERATE	-	-1	-	TACR2	HGNC	HGNC:11527	YES	ENSP00000362403	P21452	-	UPI0000061EE3	-	3/5	-	Gene3D:1.20.1070.10,Pfam_domain:PF00001,PROSITE_profiles:PS50262,hmmpanther:PTHR43919,hmmpanther:PTHR43919:SF4,SMART_domains:SM01381,Superfamily_domains:SSF81321,Conserved_Domains:cd16004	1	1	1	1	1	0.9999	0.9999	0.9999	1	0.9999	1	0.9999	1	1	1	gnomAD_ASJ,gnomAD_FIN,gnomAD_OTH,gnomAD_SAS,AFR,AMR,EAS,EUR,SAS	-	-	-	-	-	-	-

Jumbling through that, you can see the allele frequencies for `gnomAD_AF` is **0.9999**. This seems odd to me. How could this variant be a COSMIC (cancer database) missense variant, with `MODERATE` consequence, and have **99.99%** frequency. I'm lost on how to interpret this.

Maybe I am misunderstanding how gnomAD scores allele frequencies, hence posting this question here.

**Does anyone know how gnomAD allele frequencies (as outputted by Ensembl's VEP) should be interpreted?**

Answer: The reference allele, in this case, appears to be a very rare allele. The reference allele is whatever was found in the reference genome, which is a genomic region of a real person, which means that it can be a rare or private allele. In this case, the reference, A, is very rare.

The VEP is giving you the allele frequency of the alternative allele for the variant at this locus, which is [rs55953810](http://www.ensembl.org/Homo_sapiens/Variation/Population?db=core;r=10:69408429-69409429;tl=2uCd0VOAH86sXHFW-4871232;v=rs55953810;vdb=variation;vf=73998235). The alternative allele given in your VEP input is T, so the allele frequency it gives you is for T.

VEP reads VCF format as standard, assuming that the alleles are the forward strand alleles. Your specification of strand in the INFO column is ignored because this is not the standard way to write VCF. This is why your alleles have not been converted.

It appears that you have run your VEP input against GRCh38. gnomAD coordinates are on GRCh37, and are remapped onto GRCh38, so can be looked up by the VEP to find the allele frequencies. This is why there is no variant at that locus in gnomAD, but [that variant identifier and its frequencies do exist](http://gnomad.broadinstitute.org/variant/10-71168685-A-T).
Question: Suppose you have three experimental conditions for your RNA-seq: control, compound A treatment, compound B treatment, each with replicates. Your RNA-seq samples look like the following:

c1,c2,a1,a2,b1,b2

And your factors are (1,1,2,2,3,3).

And you use edgeR for the analysis. You have two questions to answer: what are the DEGs for treatment A, and what are the DEGs for treatment B. Let's just focus on the 1st question -- what are the DEGs for treatment A?

In terms of sample inclusion for analysis, we have two options:

1) c1,c2, a1,a2. You just compare (a1,a2) with (c1,c2), e.g., coef=2

2) c1,c2, a1,a2, b1,b2. You compare (a1,a2) with (c1,c2), but you use contrast=c(-1,1,0).

The pros of using option 1) is that it's simple, because we do it all the time -- pairwise comparison. The pros of option 2) is that it uses more samples and should give you more statistic power, the cons is that it's not as simple, and the inclusion of condition B samples could complicate your analysis and usually makes you hard to explain to your biologists.

I compared the above two options, and found they generated somewhat different results. Power-wise (min p-values) sometimes option 1) is better, sometimes 2) is better.

Which way is a better practice?

Thanks!

Answer: In my opinion, the second option is by far the best, for two reasons :

1. As you said, the estimation of variance is better with more samples. By the way, this doesn't always mean more power (min *pvalues*). In some cases, the added samples increase the estimation of variance, so the *pvalue* increases, but for good reason.

2. It makes sense to normalize and process all your samples together, since they are from the same study. The fact that you can answer both your questions (effect of A and B) from the same dataset, only changing the contrast, makes **the option B simpler than A**. If you choose option A, the size factors and the estimation of variance would not be the same in dataset (c1,c2, a1,a2) and (c1,c2, b1,b2). Imagine showing biologist two excel sheets with the normalized raw counts, one for dataset (A vs C), the other for dataset (B vs C). Wouldn't it be confusing for him to see that the value for gene X in condition c1 is different in both sheets ?
Question: Hi everyone!

I am working with RNA-seq paired-end data. When I run tophat-fusion on a dataset, I get mapping percentage of 70-80% for both left and right reads and the concordant pair alignment rate is about 72%. However, if I run Tophat without fusion search, I get the same mapping percentage for left and right reads, but the concodant alignment rate decreases to 30-40%. I am running both commands with the same arguments.

Does anyone know why the concordant pair alignment rate decreases with tophat-fusion? My guess would be that it should increase the concordant alignment rate as it is going to be able to map fusion events that would be unmapped by tophat.

Thank you very much!

Answer: Tophat fusion removes reads that map to more than two places in the genome, thus, the tophat mapping percentage can be higher because it won't remove multimappers by default
Question: I am building a large database for my analysis tool ~3 Terabytes.

The issue is distributing this database for people to use it.

A solution we have learned from working with sequencing cores is to ship a hard drive by FedEx. A possibly expensive and annoying proposition.

Another solution is to host the files for download, but this seems possibly expensive or possibly impractical for the audience. (On this point of practicality -- most Universities have gigabit level bandwidth, perhaps 3 Terabytes is no longer a ludicrous download size, simply a big one.) At 1 gigabit, I estimate the download would take 400 minutes.

Since this is a niche' problem, I am hoping someone may have had a similar issue and can provide advice?

Thank you
Jeremy Cox

Answer: Shipping disk is the most reliable solution once you've reached the TB range. The problem with sharing over the network is that the speed is only going to be as high as the network bottleneck between the two ends (including bottleneck inside the institutional networks) then there's some overhead depending on the protocol. Also the connections are often not stable over the many hours (up to days) required to transfer a large data set so you need mechanisms that can automatically reconnect and resume where the transfer was interrupted eventually replacing files that were corrupted when the connection was lost. There are technologies for large data transfer over the network (e.g. Aspera, GridFTP) but you end up running into costs and deployment issues.
You should also consider not transferring the data. Maybe you can provide local compute to your collaborators or allow them to extract relevant subsets and only transfer those.
Question: Hi guys,

I wonder if there're any up-to date publicly available genomes of bacteria, I could use as an input for Blast/Blat (fasta/2bit format) search. If so, would you send me links? I'm aware of database '16S ribosomal RNA sequences (Bacteria and Archaea)', however it's coverage is rather low. I appreciate your experiences.

EDIT: thank you guys, you've helped me a lot. It's a shame, that as an answer just one can be picked :)

Answer: You can get a lot from Human Microbiome Project.
Question: I am trying to subset a fasta file like described in another post (https://www.biostars.org/p/93056/). I simply want to be able to extract a subset of sequences from a fasta file.

I load the package and the data as described in the other thread:

``
library("seqinr")

fastafile<- read.fasta(file = "proteins.fasta",
seqtype = "AA",as.string = TRUE, set.attributes = FALSE)
```

as well as the list of IDs I was to use for subsetting (the column in `subsetlist` is called `id`).

subsetlist<-read.table("~/scripts/subsetfasta/test.txt", header=TRUE)

when I attempt to use the solution from the previous thread:

fastafile[names(fastafile) %in% subsetlist$id]

I get the following:

named list()

What am I doing wrong or missing?

best regards

Henrik

Answer: Looks like none of your fasta record names are in the subset list. You can check by just calling

names(fastafile) %in% subsetlist$id

Which will be a vector of FALSEs if that's the case.

You probably want to check exactly how the subset list and fasta file IDs are formatted, and use (g)sub, paste or related functions to get them to match.
Question: Hi,

I'm working on RNA-seq data and have aligned my reads to the assembled transcript using bwa-mem. When I tried to run Picard Tools (`SamFormatConverter.jar`), the SAM dictionary threw an exception

"Cannot add sequence that already exists in SAMSequenceDictionary".

I then checked for duplication in the SAM file header and sure enough, two entries were duplicated. How do I correct this? Commands used are given below:

bwa index contigs.fa
bwa mem contigs.fa reads_1.fq reads_2.fq >samFile.sam
java -Xmx4g -jar SamFormatConverter.jar INPUT=samFile.sam OUTPUT=bamFile.bam

This is where I faced the error.

I checked for duplicates using:

samtools view -HS samFile.sam | sort | uniq -c

2 @SQ entries were present >1 times.

Any help would be appreciated, thanks!

Answer: The `@SQ` headers that bwa emits come directly from the reference sequences in the index you're mapping against.

So the likely explanation for your problem is that the same entries are duplicated in your `contigs.fa` file. You can check by looking at the `>` headers in your fasta file:

grep '>' contigs.fa | sort | uniq -c

(If there are header lines like `>foo blah` and `>foo blah blah` then these will lead to duplicate `@SQ SN:foo` SAM headers while not being shown as exact duplicate lines by that `uniq` command. So you may have to be a little careful in interpreting the results of that command.)
Question: Hello dears,

We can index the files to access random sequences fast, and it should be sorted before indexed for the SAM/BAM.

But when it turns to FASTA files, I could not understand why we can directly index without sort first ?

samtools faidx

Are there any easy explanation? THANKS

Answer: Hello https://www.biostars.org/u/45640/ ,

indexing is a very fascinating topic.

The index file produced by `samtools faidx` and `.bai` have very different structures. I guess there are mainly two reasons for it:

1. In a `bam`file we have typically much more entrys than in a `fasta` file.
2. The way we query the data. For a `fasta` file we typically ask "Give me the sequence with the id XY". For `bam` files we ask "Give me all reads that overlap a region"

The `fasta` index is quite simply. It just contains the name of sequences, where in our file the header starts, how long the header is and how much bases the sequence have. See the [specs][1] for it. As the number of sequence in a `fasta` file is quite small (compared to a `bam` file) we can just iterate over the index file to find the offset of a sequence we like to have in a reasonable time.

In case of the `bam` the index file is organized in bins, which contains the offset of reads that overlap a region. See the [sam specs][2] for it . To be able to say where a bin begins and end it is necessary to sort the `bam` file.

fin swimmer

[1]: http://www.htslib.org/doc/faidx.html
[2]: http://samtools.github.io/hts-specs/SAMv1.pdf