Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Reference file for RNA Seq Analysis

    Hello,

    I have .sam files that I would like to feed into htseq-count. However, I must first have a reference file in .gtf format, I believe.

    I have access to an "assembled transcriptome" that is .fa.gz format of the same species from the .sam files. Could this possibly be used for the htseq-count as reference sequence input, and if so, should I do something to change its format as an acceptable input?

    I am sorry for the noobie question!

  • #2
    Did you just map to the assembled transcriptome? If so, you don't need to use htseq-count and can instead just use "samtools idxstats". You might need to filter by MAPQ or something to get rid of multimappers, but that'll be simpler (you'll just output this to a new file and index that).

    If you didn't map to the assemble transcriptome but instead mapped to the genome, then we'll need some more details about your experiment (in particular, the species used). htseq-count needs an annotation file because it was designed to be used for experiments where reads are mapped to the genome.

    Comment


    • #3
      Thanks you!

      Hello, I am not sure what was done to the file (I downloaded it). Is there some way I could determine that? The first 50 lines are in the format:

      >adi_EST_assem_1
      AGTAACACTGTGAAGTTGTAAAAATTCTGCATGAAAAAGCAGTGTTGTGACTCATAGTAT
      GCAAAGAACTTTCGCGAGAATCTTTTCCACAAGTTGTTGTGACTAGATTTGAAAGAGAAG
      AATGGGTATGACGTCAGTTTCAGGCGATCTACCTGATATTATGTGGAAATCAAGATGTTG
      CTGGGAGTCCCAAAATGAATACAAACATAACAATAGCTTGTCATGCCTGAGTTATTCCAA
      AGCATATGCGAACGGCAAATGGAACTCGGTAAAATGAAGCTTATCTCTTGCATGTTATAC
      CTTCTACATAGCGTCGTATTATTAGATGGATTTCAGTGGTCTCATCGACAACAGGCAAGC
      CACGTTATCAGCGTTGAAGCTGGTTCAGATGTTTCTCTTCCATGCGAATACGAACTCACG
      CCACAAGAACAGCAGGAAGCTGATGTCTTTCATCTATTAACATGGACACGAGAAGAGCCT
      CATTTCAGCGATCGGTGGACCGGATTGGCGATAAATTCGACTCTAAGTGAAAGCAAAGTT
      ATTTATGATGATCCACAGCGAATTTTCATCGTTGATGGAACACTTACTGTCAAGAACATA
      ACCGTTAAAGATAGGACTCGGTATCAATGCGCATTTCAAAGCAGTTTCTTCACGACACCC
      AGCATAATTAACTTGGACGTCAAATATCCTCCAGAGATTATGCTCAAACCATTGTCCCAA
      GATATCATTGAAGGGAATGATGTCAGGCACTGTGCAATGCATCTGGAAATCCACAGCCTA
      AGATAACTTGGACGAGGCAAGGAAACAGTGATGTCCTTTCTTCATCTGAGACCCTTTTGT
      TGCGTAATTTAGTCAGCGAAGATGATGGATCAGTGTATACGTGCAGGGTCGAAAATTACC
      TTGGATTGAAGCAAGCCTCTGTTACAATAACCATGCAATATAAGCCAAAGGACACACGAG
      TTACTATTTCATCTGGTGAGGTGAAGCTAGGAGATCGTGTGGACATCTCATGCTGGGCAC
      GTGCAAACCCTTCACCAGAGTACAAATTTTATTATAACAACAAACTGATCAGGTGGTCAA
      GTAAAGGACTGCTTTCATTAACGAGTGTTAATACTGAAGATCAAGGCACGTATCAATGTG
      TTCCAGTGAATTACCTGGGGGATGGATCTGGGGCATCTGTTACTGTAACTTTGTCCACTG
      GTGACAAGGCCTTCCCAGTATGGGCGTATGCTGCAATAGCAAGTGGTGTTCTCTTCCTTC
      TTTTAATTGCTGCGGGAATTGCATTGTGCAGGGCAAGAAAAAAATGGAAGGCAGCAGATG
      ACGGGACTAAGACAGCCCAACCGAAGGAGAGCACAATTAGCAGGCCTAGAGATGCTAAAG
      AGGCACATGTAGGTCCAGGAGATTTCTATGGAGACCACAATGGCAGTATTGAGATTCCTC
      ATCTTAATGGAGCTCTGAGTTTCAGTGACATCAGAAGTAGTCAAGAGTTCCCAGGTTTAA
      AGGGTGCTTTT
      >adi_EST_assem_2
      AGTTCAGAACGGATTGCGTGGTTTTGGCAAAAGATAGAACAGCTGAGTTGTGAAAGGAAA
      TTTGGAGTAAAGTAATTGGAATGGGTAGAAAGAACAGCAGATGTAGCCTAGCTGGGGTCG
      TCAGCCTGAGCGTACTCTTCTGCATCGTAAGCTTTTTCATGATCGTTGACGCATCACATC
      AAAGTAAACGAACACAGAGGCGAGAGAGTCTCGACAATAACAGAGTGAACTTGAGCGCCG
      GAAATTTCCTCGTAGAATGGAACGTTGACAAAAACTCTGGAAGCATTGAATTCATATTGA
      AGGTTACAGTTGATGAAAAAGGCTGGATTTTCTTTGGATTCATGCCCATCTCAATTAACA
      GTTCCACGAAGCCAAGCGGCATTGAAGATTCTACAGGAGACTTCGTGGTAACATGGTTAT
      TACCTGGCTCTAACGACAGAGAAACTTTGGATTTGAATACAATAAAAAAGAAAGGGGAGC
      TAGGGATAGACACCGACAACAATTATGGGATTACAGCAGGGGATTCTAACAAGAATTCAA
      GTGTACAGGTGCTTCACATTACCAGGAAACTTGACACAGGCGACCCACAAGATGTGCCCA
      TCACAAACCAACCAATGTGGATGTTTGCTGCATGGGGAACATCTGAAGAATTCACAAAAA
      ACTTCACCACAATTAAGGACAACATAACCAACTCCCTAGTGGTGGTAAAAGTGGTTTTTT

      I'm guessing that does not help, though, right?

      Where I downloaded the file, it said "Assembled Transcriptome v1.0, (compressed with gzip), total 36,780 contigs for 29,364,984 bp"

      I am doing this just for practice to get a feel for RNASeq. But I have limited time, so I am not too concerned about quality at all (NOT to be published etc)!!

      Comment


      • #4
        Clarification

        dpryan:

        Thanks again.

        I realize I did not show what the .sam files look like. The head is as such:

        @HD VN:1.0 SO:coordinate
        @SQ SN:scaf1 LN:42540
        @SQ SN:scaf10 LN:28213
        @SQ SN:scaf100 LN:68029
        @SQ SN:scaf1000 LN:50671
        @SQ SN:scaf10000 LN:2123

        and tail is as such:

        WUSI-EAS1631:100:FC63B34AAXX:1:99:14443:9257 0 scaf9999 4971 50 21M406N79M * 0 0 CAAAACGCTGACGCCATCCACCTCAACCAAATGGTCTCCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGC HHHHHHFHDHHHHHHHEHHHHHBHGHHHHHHHFHEDFFGFHBHFHCFED>BBDDBCEC@B>@>CADA9B@ECB@BEAEBA<@@@@############### XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1
        HWUSI-EAS1631:100:FC63B34AAXX:8:13:5760:14033 16 scaf9999 5406 50 87M508N13M * 0 0 AATGGTCTCCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTG GGEDGCEBEEHGHEBIIEIHDIFHIIGHIHHAFHIEIIIIGIIIFIHHIHIHIIIIIHIIIIIGIIIIIIIIIIII@IIIIIIIIIHIIIIIIIIIIIII XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1
        HWUSI-EAS1631:100:FC63B34AAXX:2:68:12267:10621 0 scaf9999 5414 50 79M508N21M * 0 0 CCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCCTCTAATGAATTTTGAGTTCCTT HHHHHHHHHHHHHHDHHHHDDGGDGHHGHFHFFFHGDGDFHHHHHGH@CGEBHFEGDHEBBBDEC>B7CBBBB@B64*9==9>6==@=BBB==<===4;= XA:i:1 MD:Z:77A22 NM:i:1 XS:A:- NH:i:1
        HWUSI-EAS1631:100:FC63B34AAXX:2:67:17937:17163 16 scaf9999 5424 50 69M508N31M * 0 0 CGACTCCATTCTGCCCCACCGGCCCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTGAGTTCCTTCCAAACTAAT ####################@?A><EC3GD<DI@IGFIIBDHIGFGGBDEGEB@GE<EDGDBG@GGGDEGDEHFDGIGE@BDGGBDGHGFIIIGG@II XA:i:2 MD:Z:13A8A77 NM:i:2 XS:A:- NH:i:1
        HWUSI-EAS1631:100:FC63B34AAXX:1:99:14443:9257 16 scaf9999 5431 50 62M508N38M * 0 0 ATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTGAGTTCCTTCCAAACTAATTCCTAAG >B>?8G5DBFC<BFDFBIDE>BG@BGB>BGGEDGCEE?C?EAGDHEDHIIIIHIGBGGGGDGG8GGGBGHDIIHGIFFDGGBGGBGGGGGGGGBGEGEGE XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1

        Could that indicate whether the .sam files have been mapped? If so, then are you saying I can run "samtools idxstats" on the .sam files?

        I did try to run samtools idxstats on them using the command:

        ../../samtools-0.1.7_i386-darwin/samtools idxstats adu.sam

        But get the error:

        [main] unrecognized command 'idxstats'

        However, samtools did seem to be working when I earlier converted between files using samtools views, as such:
        ../samtools-0.1.7_i386-darwin/samtools view -h -o adu.sam adu.bam

        I'm not sure what might be causing that error (yes, I am very new to this)...sorry (

        Thanks so much again, for any help......!!!

        Comment


        • #5
          One other detail

          And the species is coral (A. digitifera)

          Comment


          • #6
            You can probably use the gff file. The genome sequence associated with it is here, in case it's different from the one you aligned against.

            Comment


            • #7
              Thank you!

              dpryan:

              Thank you so much!

              I was wondering, where did you find this .gff file?

              I followed steps to successfully install htseq-count. I believe it is actually working (rare for me), as I can type into the terminal htseq-count, and it auto-fills as a command.

              However, I tried a simple command like:

              htseq-count adu.sam aug_repeatmask_pasa_input.gff3 > ADUcount.out

              and got the error:

              Feature aug_v2a.20244.t1.exon1 does not contain a 'gene_id' attribute
              [Exception type: ValueError, raised in count.py:53]

              I am assuming the .gff3 file is missing a "gene_id" field, but am not sure how to fix that?

              Sorry for another question, and thanks for your help with this!

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Strategies for Sequencing Challenging Samples
                by seqadmin


                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                03-22-2024, 06:39 AM
              • seqadmin
                Techniques and Challenges in Conservation Genomics
                by seqadmin



                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                Avian Conservation
                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                03-08-2024, 10:41 AM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 06:37 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, Yesterday, 06:07 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-22-2024, 10:03 AM
              0 responses
              51 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-21-2024, 07:32 AM
              0 responses
              67 views
              0 likes
              Last Post seqadmin  
              Working...
              X