gpertea
STAR -> StringTie run giving very small number of novel exons/introns/loci (<1.0% of total)
I am running STAR and StringTie on RNAseq data from 4 samples (PE 100bp, 60M reads/sample, stranded) and am looking at the merged.stats file. The numbers for novel stuff being found in my samples is REALLY low, which makes me think I need to change my run parameters somewhere to be less stringent somehow. Here is all the code I ran and the contents of my merged stats file:
#Generate genome index using STAR STAR --runThreadN 16 --runMode genomeGenerate --genomeDir star --genomeFastaFiles genome/GRCm38.primary_assembly.genome.fa --sjdbOverhang 99 --sjdbGTFfile genes/gencode.vM25.primary_assembly.annotation.gtf
#Align to genome using STAR STAR --genomeDir star_THIS --readFilesCommand zcat --readFilesIn samples/2_Forward.fq.gz samples/2_Reverse.fq.gz --outSAMtype BAM SortedByCoordinate --outFilterType BySJout --alignSJoverhangMin 8 --limitBAMsortRAM 16000000000 --outSAMunmapped Within --twopassMode Basic --quantMode TranscriptomeSAM GeneCounts --outSAMstrandField intronMotif --runThreadN 16 --outFileNamePrefix "2_star_THIS/" STAR --genomeDir star_THIS --readFilesCommand zcat --readFilesIn samples/3_Forward.fq.gz samples/3_Reverse.fq.gz --outSAMtype BAM SortedByCoordinate --outFilterType BySJout --alignSJoverhangMin 8 --limitBAMsortRAM 16000000000 --outSAMunmapped Within --twopassMode Basic --quantMode TranscriptomeSAM GeneCounts --outSAMstrandField intronMotif --runThreadN 16 --outFileNamePrefix "3_star_THIS/" STAR --genomeDir star_THIS --readFilesCommand zcat --readFilesIn samples/4_Forward.fq.gz samples/4_Reverse.fq.gz --outSAMtype BAM SortedByCoordinate --outFilterType BySJout --alignSJoverhangMin 8 --limitBAMsortRAM 16000000000 --outSAMunmapped Within --twopassMode Basic --quantMode TranscriptomeSAM GeneCounts --outSAMstrandField intronMotif --runThreadN 16 --outFileNamePrefix "4_star_THIS/" STAR --genomeDir star_THIS --readFilesCommand zcat --readFilesIn samples/5_Forward.fq.gz samples/5_Reverse.fq.gz --outSAMtype BAM SortedByCoordinate --outFilterType BySJout --alignSJoverhangMin 8 --limitBAMsortRAM 16000000000 --outSAMunmapped Within --twopassMode Basic --quantMode TranscriptomeSAM GeneCounts --outSAMstrandField intronMotif --runThreadN 16 --outFileNamePrefix "5_star_THIS/"
#Assemble Samples 2 to 5 only stringtie-2.1.4/stringtie -p 40 -G genes/gencode.vM25.primary_assembly.annotation.gtf -o WT2_star_THIS.gtf 2_star_THIS/Aligned.sortedByCoord.out.bam stringtie-2.1.4/stringtie -p 40 -G genes/gencode.vM25.primary_assembly.annotation.gtf -o WT3_star_THIS.gtf 3_star_THIS/Aligned.sortedByCoord.out.bam stringtie-2.1.4/stringtie -p 40 -G genes/gencode.vM25.primary_assembly.annotation.gtf -o WT4_star_THIS.gtf 4_star_THIS/Aligned.sortedByCoord.out.bam stringtie-2.1.4/stringtie -p 40 -G genes/gencode.vM25.primary_assembly.annotation.gtf -o WT5_star_THIS.gtf 5_star_THIS/Aligned.sortedByCoord.out.bam
#Merge sample 2 to 5 only with main annotation stringtie-2.1.4/stringtie --merge -p 40 -G genes/gencode.vM25.primary_assembly.annotation.gtf -o adipo_stringtie_merged_2to5only.gtf mergelist_2to5only.txt Compare merged vs main annotation gffcompare-0.12.1/gffcompare -r genes/gencode.vM25.primary_assembly.annotation.gtf -G -o merged_2to5only adipo_stringtie_merged_2to5only.gtf #Contents of merged_2to5only.stats file
gffcompare v0.12.1 | Command line was: gffcompare-0.12.1/gffcompare -r genes/gencode.vM25.primary_assembly.annotation.gtf -G -o merged_2to5only adipo_stringtie_merged_2to5only.gtf Summary for dataset: adipo_stringtie_merged_2to5only.gtf
Query mRNAs : 148525 in 53729 loci (121240 multi-exon transcripts) (20486 multi-transcript loci, ~2.8 transcripts per locus) Reference mRNAs : 142004 in 53534 loci (115073 multi-exon) Super-loci w/ reference transcripts: 53294 -----------------| Sensitivity | Precision | Base level: 100.0 | 96.6 | Exon level: 94.1 | 97.3 | Intron level: 100.0 | 98.9 | Intron chain level: 99.9 | 94.8 | Transcript level: 99.4 | 95.0 | Locus level: 99.7 | 99.0 |
Matching intron chains: 114995 Matching transcripts: 141172 Matching loci: 53363
Missed exons: 0/446814 ( 0.0%)
Novel exons: 1754/421810 ( 0.4%)
Missed introns: 4/284950 ( 0.0%)
Novel introns: 783/288229 ( 0.3%)
Missed loci: 0/53534 ( 0.0%)
Novel loci: 435/53729 ( 0.8%)
Total union super-loci across all input datasets: 53729 148525 out of 148525 consensus transcripts written in merged_2to5only.annotated.gtf (0 discarded as redundant)
I also noticed in some recent simulated experiments that, with the default parameters, aligning to a genome index without any junction data, STAR has a much lower sensitivity/recovery rate for the spliced alignments (around 74%) compared to HISAT2 (about 97%).
When junction data is present STAR is going to find those alignments/isoforms of course -- but it has a much lower chance of finding any "novel" junctions (and thus novel isoforms), compared to HISAT2.
This does not have much to do with StringTie itself though. Switching to HISAT2 as your aligner (without the --dta option!) would increase the chances of novel isoform discovery with StringTie.
