This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Low alignment with Trinity but high alignment in SAM files

Hello,

I have performed alignment with Trinity and got only around 2% alignment. My SAM files and trimmed reads seemed to be fine; the SAM files had high alignment rate based on what was returned to me with samtools. Is this possible or am I doing something wrong with Trinity?

I appreciate the help

alignment genome rna-seq next-gen trinity
• 4,611 views • link
updated by h.mon h.mon • written by Moneeb Bajwa • Moneeb Bajwa

Could you please share the commands used for both alignment as well for getting stat from SAM.

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

I just used gsnap for the SAM files by aligning trimmed mouse reads to the mouse genome (GRCm38.primary_assembly.genome.fa from GENCODE). I used gmap to build the database first and built intron index with the gencode.vM18.primary_assembly.annotation.gff3 file. I used this for the SAM stat: 

#Count the total number of reads in the alignment
total="$(samtools view -c $samFile)"
#Count the aligned reads in the alignment
mapped="$(samtools view -F4 -c $samFile)"

For Trinity I did this (they were single-end, strand-specific forward reads):

nice -n19 $trinityPath/util/align_and_estimate_abundance.pl \
--transcripts Trinity-GG.fasta --aln_method bowtie2 --est_method eXpress \
--trinity_mode --output_dir $outPath$sample --seqType fq --SS_lib_type F \
--single $Reads --thread_count 1

And for the Trinity-GG.fasta I did this: 

Trinity --genome_guided_bam MergedBams.bam \
--genome_guided_max_intron 10000 \
--max_memory 10G --CPU 4
• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

What is your library prep kit? Is it really stranded? Is it from Illumina? In general, stranded Illumina kits use dUTP method, with Trinity, you should use --SS_lib_type R for single-end reads, and --SS_lib_type RF for paired-end reads.

Are the reads you are trying to align the same used to assemble the transcriptome?

Are you working with Mus musculus? Do you really need to assemble a transcriptome?

• 0 views • link
written by h.mon h.mon

Here are the reads' link which gives all the info needed: https://www.ncbi.nlm.nih.gov/bioproject/PRJDB5631. I need to do differential expression analysis so I need to assemble the transcriptome right? Also the merged BAM file was only like 10G while the SAM files in total were around 150G; is that normal?

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

Yea actually it is still not working after rerunning the sorting and merging; the BAM files are quite small still and I believe that is messing everything up

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

1 answer

First of all I will answer the question you asked on a comment, but it is the most important question here:

I need to do differential expression analysis so I need to assemble the transcriptome right?

You do not need to assemble the transcriptome to perform gene / transcript differential expression analysis. There is a very, very good genome assembly for mouse. Also, there is a very, very good annotation for the mouse genome - and consequently, there is a very good transcriptome for the mouse. You do not need to assemble, it is better to map / count using the existing reference genome / transcrptome. In addition, the reads from the project you linked are single-end and only 36bp long, they will not provide a good transcriptome assembly.

You can perform gene-level differential expression by first mapping the reads to the genome with STAR to get gene expression counts (as I have told you before here), and then using DESeq2 or edgeR to perform the differential expression statistical analysis.

You can perform transcript-level differential expression quantifying transcript counts with Salmon or kallisto, and performing the statistical analysis with sleuth.

Now, for your original question:

I did map one of the files to the mouse genome with STAR, obtaining gene counts:

Sox17   9368    81      9287
Mrpl15  4473    30      4443
Lypla1  3122    20      3102
Tcea1   397     6       391

One can use these counts to infer strandedness of the reads. I know the protocol is stranded, so the first column of counts is automatically excluded, as it represents unstranded counts (check STAR manual). It is easy to see then that the true counts are in the fourth column, which the STAR manual describes as:

column 4:  counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)

This means the library prep kit produces what we call "reverse" stranded reads. The Trinity manual says the correct parameter to indicate this library type is --SS_lib_type R, which is the opposite of what you are using in your call to the align_and_estimate_abundance.pl script. So maybe if you use align_and_estimate_abundance.pl --SS_lib_type R you will get good mapping rate.

But remember, you do not need to assemble the transcriptome.

• 0 views • link
written by h.mon h.mon

I also did try --SS_lib_type R, but it did not change anything for Trinity. I will just go ahead with the STAR option you mentioned. But I still am not sure why the STAR method gives reverse when the link says forward clearly. Maybe they are talking about two different things?

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

Is there another route I can go without building the transcriptome since STAR is not working?

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

Hi,

STAR alignment happens in two steps the first time and in one step the subsequent times. Like Tophat, Bowtie etc has index created prior to aligment, so does STAR, but here instead of index file, a whole directory with multiple files will be created whose name will be taken from genome_dir argument. Once it is created, you may use it like you have used. Also please read about outSAMtype. I dont think None is the best you can give there.

In short,

In Step1use args --runMode genomeGenerate providing fasta and gtf alone(fastqs do not come into picture here). Once that is done, go for the other runmode alignReads giving your fastqs and genome dir created in the previous step. STAR is not very complicated. Don't worry. Try this and for specifics go through the STAR manual. It's very helpful.

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

I get this error now: EXITING because of FATAL PARAMETER ERROR: limitGenomeGenerateRAM=31000000000 is too small for your genome. SOLUTION: please specify --limitGenomeGenerateRAM not less than 95444909440 and make that much RAM available.

Isn't that 95TB? I can't do that.

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

This is what I did:

#!/bin/sh
STAR --runThreadN 5 --runMode genomeGenerate --genomeDir genome_dir \
--genomeFastaFiles gencode.vM18.transcripts.fa \
--sjdbGTFfile gencode.vM18.primary_assembly.annotation.gtf \
1>prepStar.log 2>prepStar.err &
• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

Ok but I keep getting this message: Fatal INPUT FILE error, no valid exon lines in the GTF file: gencode.vM18.annotation.gtf. Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA file.

I am using the GTF and fasta files from GENCODE so I don't understand why there is an issue. I tried using the two labelled PRI and I also tried using the first one of each in this link: https://www.gencodegenes.org/mouse_releases/current.html

• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

This is the code I used:

#!/bin/sh
STAR --runThreadN 5 --runMode genomeGenerate --genomeDir genome_dir \
--genomeFastaFiles gencode.vM18.transcripts.fa \
--sjdbGTFfile gencode.vM18.annotation.gtf \
--limitGenomeGenerateRAM 95444909440 \
1>prepStar.log 2>prepStar.err &
• 0 views • link
written by Moneeb Bajwa • Moneeb Bajwa

Log in to answer this question.

More posts like this