Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Reads are 0 in htseq-count

    Hello,

    I have a trouble using htseq-count, I used TopHat to aling my data, but when I tried to count the number of reads all the genes present 0 read.

    I checked the chromosome names in ./gtf file: cut -f 1 -d " " ~/HG19/genes.gtf | cut -f 1 | sort | uniq So all of them are started with "chr". The same result was accepted by bowtie2-inspect --names ~/HG19/genome. (in .gtf file, I've seen the numbers like "chrHSCHR9_2_CTG3" and don't know what it is).

    1) programs/tophat2 -p 4 -r 50 -G ~/HG19/genes.gtf -o tst_output HG19/genome tst_fastq Result: accepted_hits.bam deletions.bed junctions.bed prep_reads.info align_summary.txt insertions.bed logs unmapped.bam

    samtools view -H ./accepted_hits.bam:
    @HD VN:1.0 SO:coordinate @SQ SN:chr1 LN:249250621 @SQ SN:chr10 LN:135534747 @SQ SN:chr11 LN:135006516 @SQ SN:chr12 LN:133851895 @SQ SN:chr13 LN:115169878 @SQ SN:chr14 LN:107349540 @SQ SN:chr15 LN:102531392 @SQ SN:chr16 LN:90354753 @SQ SN:chr17 LN:81195210 @SQ SN:chr18 LN:78077248 @SQ SN:chr19 LN:59128983 @SQ SN:chr2 LN:243199373 @SQ SN:chr20 LN:63025520 @SQ SN:chr21 LN:48129895 @SQ SN:chr22 LN:51304566 @SQ SN:chr3 LN:198022430 @SQ SN:chr4 LN:191154276 @SQ SN:chr5 LN:180915260 @SQ SN:chr6 LN:171115067 @SQ SN:chr7 LN:159138663 @SQ SN:chr8 LN:146364022 @SQ SN:chr9 LN:141213431 @SQ SN:chrM LN:16571 @SQ SN:chrX LN:155270560 @SQ SN:chrY LN:59373566 @PG ID:TopHat VN:2.1.1 CLrograms/tophat -p 4 -r 50 -G ~/HG19/genes.gtf -o tst_output HG19/genome tst.fastq
    less align_summary.txt:
    Reads: Input : 99 Mapped : 76 (76.8% of input) of these: 1 ( 1.3%) have multiple alignments (0 have >20) 76.8% overall read mapping rate.
    2) samtools view -h tst_output/accepted_hits.bam | htseq-count - ~/HG19/genes.gtf > tst_count.txt tail tst_count.txt:
    ENST00000610276 0 ENST00000610277 0 ENST00000610278 0 ENST00000610279 0 ENST00000610280 0 __no_feature 40 __ambiguous 27 __too_low_aQual 0 __not_aligned 0 __alignment_not_unique 15
    I have no idea what should I do, so could you please help me.

    PS: after samtools sort -n -o temp-sorted.bam accepted_hits.bam I have the same result ("SO:queryname" in temp-sorted.bam file).

    Thanks,
    Maria.

  • #2
    It looks like you are mixing and matching genome builds and annotation. "chr" designations on your chromosomes indicates that the data is from UCSC but your annotation has ENST* identifiers which indicates Ensembl as the origin.

    I suggest that you find annotation that uses the chr designations for doing the counting.

    BTW: Please stop using TopHat. It is not even recommended by its authors. There are current options like STAR/bbmap/HISAT2 that should be used instead.

    Comment


    • #3
      Originally posted by GenoMax View Post
      It looks like you are mixing and matching genome builds and annotation. "chr" designations on your chromosomes indicates that the data is from UCSC but your annotation has ENST* identifiers which indicates Ensembl as the origin.

      I suggest that you find annotation that uses the chr designations for doing the counting.

      BTW: Please stop using TopHat. It is not even recommended by its authors. There are current options like STAR/bbmap/HISAT2 that should be used instead.
      Thanks, GenoMax, you are absolutely right: my annotation was Ensembl and my .fa was UCSC.

      I downloaded the annotation from http://genome.ucsc.edu/cgi-bin/hgTables - is it correct?

      head my .gtf file:
      chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
      chr1 hg19_knownGene exon 12613 12721 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
      chr1 hg19_knownGene exon 13221 14409 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
      All reads are zero in htseq-count again.

      my .fastq file:
      @ERR599187.1 HWI-ST1034:1070DURACXX:4:1101:1451:1992/1
      ATCGCCTCCTGNTAGGCCAGCTCTACATCTTCTTGGGTGACCACTAGTTTAATCTTATTCCATTTCTCAGGGATATCCAGCCACTGGTGCACAGCCACTGC
      +
      @@@DDDDDHHH#2AEHG>EG?GAGA>FECHIIIICEG?GFBBCGHDFIFGIIIIGGHAFHIIIGG=DDAHBHH@DEDEECBCC@BCCCCCCCBBBCCCC
      @ERR599187.2 HWI-ST1034:1070DURACXX:4:1101:1687:1994/1
      AAAAGAAGTAANATAAGAAGGCAATGCTTGTGGAATGTACAGTGCATATTGGCGGCGCACGCCTCATTACGATTCGCCTGCTTTCTTCTCCTGTTCAATCG
      +
      @@@DDDDDHHH#3AFHGGGGII@EHIGIIIIIGIBGHHIFGI0BFBGIII@H<FE8=/398=89@@CD#################################
      @ERR599187.3 HWI-ST1034:1070DURACXX:4:1101:1636:1996/1
      CACCCACGGAANCGAGAAAGAGCTATCAATCTGTCAATCCTGTCCGTGTCCGGGCCGGGTGAGGGTTCCCGTGTTGGGTCAAATTAAGCCGCAGGCTCCAC
      +
      @@@?DD=DD<F#222CFHIGIC;FHEDDGH@H4?<9BGHE<FFHGGIGIIGHA4=>C@#########################################
      tail test_count.txt
      uc031tkv.1 0
      uc031tkw.1 0
      uc031tkx.1 0
      uc031tky.1 0
      uc031tkz.1 0
      __no_feature 4
      __ambiguous 0
      __too_low_aQual 0
      __not_aligned 0
      __alignment_not_unique 0
      BTW: I see what you mean, I know about STAR and HISAT-2. It's just my task.

      Comment


      • #4
        Can you post a few lines from your SAM/BAM file? Just part of the SAM header should be fine.

        Comment


        • #5
          Thanks for your attention to my problem, GenoMax.

          samtools view -H temp.samheader.sam
          @HD VN:1.0 SO:unsorted
          @SQ SN:0_2.65543%_cov_993_len_552 LN:552
          @SQ SN:1_2.28932%_cov_1093_len_446 LN:446
          samtools view -H test_genome.bwt.samheader.sam
          @HD VN:1.0 SO:coordinate
          @SQ SN:0_2.65543%_cov_993_len_552 LN:552
          @SQ SN:10_0.802277%_cov_294_len_562 LN:562
          samtools view -H accepted_hits.bam
          @HD VN:1.0 SO:coordinate
          @SQ SN:chr1 LN:249250621
          @SQ SN:chr10 LN:135534747
          I guess I downloaded the wrong annotation because I see this in the sample-file sample_count.txt:
          A1BG 226
          A1BG-AS1 202
          A1CF 5
          However, my .gtf file doesn't contain such names of genes. I can't find the annotation hg19 from UCSC which include such names.

          Comment


          • #6
            Best solution is to get the sequence/annotation bundle from Illumina's iGenomes site (note: link is for big file download). In that download you will find a genes.gtf file in Annotation/Genes folder. Use that for your counting.

            I also recommend that you use featureCounts. Much faster and easier to use.

            Comment

            Latest Articles

            Collapse

            • seqadmin
              Exploring the Dynamics of the Tumor Microenvironment
              by seqadmin




              The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
              07-08-2024, 03:19 PM
            • seqadmin
              Exploring Human Diversity Through Large-Scale Omics
              by seqadmin


              In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
              06-25-2024, 06:43 AM

            ad_right_rmr

            Collapse

            News

            Collapse

            Topics Statistics Last Post
            Started by seqadmin, 07-19-2024, 07:20 AM
            0 responses
            139 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 07-16-2024, 05:49 AM
            0 responses
            116 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 07-15-2024, 06:53 AM
            0 responses
            109 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 07-10-2024, 07:30 AM
            0 responses
            43 views
            0 likes
            Last Post seqadmin  
            Working...
            X