1
$\begingroup$

I have mapped reads to the yeast (R64-1-1) genome with STAR like this:

STAR --runThreadN 44 --runMode alignReads --genomeDir /home/user/Documents/scripts/pyseqtools/index/STAR/R64-1-1 --readFilesIn /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R1_001_val_1.fq.gz /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R2_001_val_2.fq.gz --readFilesCommand zcat --quantMode TranscriptomeSAM GeneCounts --twopassMode Basic --outSAMunmapped None --outSAMattrRGline ID:WT_3 PU:WT_3 SM:WT_3 LB:unknown PL:illumina --outSAMtype BAM Unsorted --outTmpDir /home/user/Documents/analyses/TT-Seq/user/tmp --outFileNamePrefix /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3

I wanted to count the reads using htseq-count as follows:

htseq-count --quiet --format bam --stranded yes /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3Aligned.out.bam /home/user/Documents/scripts/pyseqtools/gtf/R64-1-1/Saccharomyces_cerevisiae.R64-1-1.106.gtf > /home/user/Documents/analyses/TT-Seq/user/htseq-count/WT_3.count.txt

The problem is that most reads map to no_feature:

__no_feature    4462786
__ambiguous 12644
__too_low_aQual 0
__not_aligned   4907
__alignment_not_unique  12529444

I have check whether the chromosome names in the BAM file and the GTF file matched (which they seem to do):

BAM:
@HD VN:1.4
@SQ SN:I    LN:230218
@SQ SN:II   LN:813184
@SQ SN:III  LN:316620
@SQ SN:IV   LN:1531933
@SQ SN:V    LN:576874
@SQ SN:VI   LN:270161
@SQ SN:VII  LN:1090940
@SQ SN:VIII LN:562643
@SQ SN:IX   LN:439888
@SQ SN:X    LN:745751
@SQ SN:XI   LN:666816
@SQ SN:XII  LN:1078177
@SQ SN:XIII LN:924431
@SQ SN:XIV  LN:784333
@SQ SN:XV   LN:1091291
@SQ SN:XVI  LN:948066
@SQ SN:Mito LN:85779

GTF: 
IV  sgd gene    8683    9756    .   -   .   gene_id "YDL246C"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL246C_mRNA-E1";
IV  sgd CDS 8686    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL246C";
IV  sgd start_codon 9754    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  8683    8685    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    17577   18566   .   -   .   gene_id "YDL243C"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL243C_mRNA-E1";
IV  sgd CDS 17580   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL243C";
IV  sgd start_codon 18564   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  17577   17579   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    1248154 1249821 .   -   .   gene_id "YDR387C"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDR387C_mRNA-E1";
IV  sgd CDS 1248157 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDR387C";
IV  sgd start_codon 1249819 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  1248154 1248156 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";

Edit: the output of ReadsPerGene.out.tab is (thanks swbarnes2) :

N_unmapped  73359652    73359652    73359652
N_multimapping  15789900    15789900    15789900
N_noFeature 192619  5872967 380420
N_ambiguous 648956  13479   430527

Would anyone know why most reads map to no_feature?

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ You are already running it with GeneCounts, which should output the same as htseq-count. What does your ReadsPerGene.out.tab say? $\endgroup$
    – swbarnes2
    Jul 27, 2022 at 22:16
  • 2
    $\begingroup$ Do your gtf contains exon entris (i.e. with "exon" on the 3rd column)? Could you please update your question with a few exon entries from your gtf? $\endgroup$
    – bepoli
    Jul 29, 2022 at 11:18

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.