STAR比对SMART-seq单细胞数据，外显子/内含子比对结果一致问题求助

Hey Leon, let's break down why your intron alignment is spitting out the exact same counts as your exon one—this is a common pitfall with how STAR interprets the parameter you're using, and we can fix it with a few straightforward adjustments.

Why Your Current Setup Isn't Working #

The core issue is the --sjdbGTFfeatureExon intron flag. This parameter doesn't restrict reads to align only to introns—instead, it tells STAR to treat the intron features in your GTF as if they were exons when building its splice junction database. At the end of the day, STAR is still aligning reads across the entire genome, and when it runs GeneCounts, it's tallying reads per gene (not per feature type). That's why your exon and intron counts are identical.

Also, quick sanity check: make sure the introns you generated with construct_introns are actually labeled as intron in your GTF (not some other feature type) and are properly linked to their parent genes.

Two Reliable Solutions to Get Spliced/Unspliced Counts #

Option 1: Use Feature-Specific GTFs + STAR #

This approach creates separate indexes and alignments tailored to exons and introns, ensuring clean counts:

Step 1: Create Exon-Only and Intron-Only GTFs

First, filter your GTFs to keep only the features we care about:

• Exon-only GTF (from original Ensembl GTF):

  awk '$3 == "exon"' GRCm38build100.gtf > GRCm38build100_exons_only.gtf
  

• Intron-only GTF (from your modified intron-included GTF):

  awk '$3 == "intron"' GRCm38build100withIntrons.gtf > GRCm38build100_introns_only.gtf
  

Step 2: Build Separate STAR Indexes (Optional but Recommended)

Building dedicated indexes avoids cross-contamination from splice junctions meant for the other feature type:

• Exon index:

  STAR \
  --runMode genomeGenerate \
  --runThreadN 8 \
  --sjdbOverhang 50 \
  --genomeDir GRCm38build100_exon_index/ \
  --genomeFastaFiles GRCm38build100.fa \
  --sjdbGTFfile GRCm38build100_exons_only.gtf
  

• Intron index:

  STAR \
  --runMode genomeGenerate \
  --runThreadN 8 \
  --sjdbOverhang 50 \
  --genomeDir GRCm38build100_intron_index/ \
  --genomeFastaFiles GRCm38build100.fa \
  --sjdbGTFfile GRCm38build100_introns_only.gtf
  

Step 3: Run STAR Alignments for Each Feature

• Exon alignment (tweaked to use the dedicated index/GTF):

  STAR \
  --runMode alignReads \
  --runThreadN 8 \
  --genomeDir GRCm38build100_exon_index/ \
  --readFilesIn $R1 \
  --outSAMtype None \
  --quantMode GeneCounts \
  --sjdbGTFfile GRCm38build100_exons_only.gtf \
  --outFileNamePrefix output/${SAMPLE}_exon_
  

• Intron alignment:

  STAR \
  --runMode alignReads \
  --runThreadN 8 \
  --genomeDir GRCm38build100_intron_index/ \
  --readFilesIn $R1 \
  --outSAMtype None \
  --quantMode GeneCounts \
  --sjdbGTFfile GRCm38build100_introns_only.gtf \
  --outFileNamePrefix output/${SAMPLE}_intron_
  

Note: I removed --twopassMode Basic here—since we're targeting specific regions, we don't need STAR to discover novel splice junctions, which will speed up runs.

Option 2: Single Alignment + featureCounts (Faster) #

If you don't want to rebuild indexes, do one full genome alignment and use featureCounts to split counts by feature type:

Step 1: Run Full Genome Alignment

STAR \
--runMode alignReads \
--runThreadN 8 \
--genomeDir GRCm38build100/50bp/ \
--readFilesIn $R1 \
--outSAMtype BAM SortedByCoordinate \
--twopassMode Basic \
--outFileNamePrefix output/${SAMPLE}_full_

Step 2: Count Exon Reads with featureCounts

featureCounts \
-T 8 \
-t exon \
-g gene_id \
-a GRCm38build100.gtf \
-o output/${SAMPLE}_exon_counts.txt \
output/${SAMPLE}_full_Aligned.sortedByCoord.out.bam

Step 3: Count Intron Reads with featureCounts

featureCounts \
-T 8 \
-t intron \
-g gene_id \
-a GRCm38build100withIntrons.gtf \
-o output/${SAMPLE}_intron_counts.txt \
output/${SAMPLE}_full_Aligned.sortedByCoord.out.bam

This is my go-to for SMART-seq data since it cuts down on redundant alignment steps.

Quick Sanity Checks #

1. Verify your intron GTF is valid:

   library(rtracklayer)
   intron_gtf <- import("GRCm38build100withIntrons.gtf")
   table(intron_gtf$type) # Should show a non-zero count for "intron"
   

2. For SMART-seq (full-length), your intron counts should be much lower than exon counts—if they're still identical, double-check that your intron-only GTF doesn't accidentally include exons.

内容的提问来源于stack exchange，提问作者Leon