RNA-seq : Quantification problem after mapped with BBMap

Question

[k.kathirvel93]

Hi EveryOne,

I have used BBMap for RNA-seq reads mapping. After mapping, i tried HTSeq, FeatureCounts and STRINGtie for quantification. But, i didn't get any results from featureCounts and HTSeq (readcount for all genes 0), STRINGtie showing that " your gtf file doesn't matched with your reference genome ". The strange is the same reference genome and annotation gtf files were tried with STAR aligner, but i didn't get any errors , working fine. Any suggestions, Thanks .

Here is the command used for BBMap :

./bbmap.sh threads=16 in1=/home/1.fastq in2=/home/2.fastq out=/home/output.sam ref=/home/Homo_sapiens.GRCh38.77.fa nodisk

[Tm]

please list some header of sequences and screenshot of gtf file to give us better idea.

Answer 1

[michael.ante]

Hi k.kathirvel93,

Next to trausch's suggestion that bbmap's cigar string is per default in sam 1.4 notation (which can be switched by setting sam=1.3), BBmap uses the full fasta header for the chromosome name. Instead of having just e.g. "21" in the sam file, you have "21 dna:chromosome chromosome:GRCh38:21:1:46709983:1 REF".

To switch this behaviour of, use the parameter trimreaddescriptions=t.

Cheers,

Michael

[k.kathirvel93]

where should I use this trimreaddescriptions=t option? Thanking you

[michael.ante]

AFAIK, you can use it either during index generation or during mapping.

Answer 2

[trausch]

I think, BBMap uses the new CIGAR characters = (sequence match) and X (sequence mismatch) and that might not be supported yet in your quantification methods. The latest version of Alfred v0.1.8 can handle these and also does RNA-Seq quantification:

alfred count_rna -g Homo_sapiens.GRCh38.77.gtf.gz <input.bam>