Announcement

Collapse
No announcement yet.

htseq-count for sam and gff3

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #16
    Almost...

    Hi Simon.
    First of all, thank you very much!!!
    Second:
    Apparently, it worked. It was an ; extra in a pair of lines so it continues reading everything just until it arrives to the #fasta part:

    htseq-count -s no -i ID s_1_sequence_clipped_tophat.bam TbruceiTreu927_TriTrypDB-3.3.3.gff
    Error occured in line 46748 of file TbruceiTreu927_TriTrypDB-3.3.3.gff.
    Error: need more than 1 value to unpack
    [Exception type: ValueError, raised in __init__.py:214]

    And when I look at this line this is what I can see:


    tryp_XI-1036e06.p1k.snoRNA.0004;Dbxref=ApiDB:tryp_XI-1036e06.p1k.snoRNA.0004,taxon:185431
    apidb|tryp_XI-1036e06.p1k ApiDB exon 22123 22200 . + . ID=apidb|exon_tryp_XI-1036e06.p1k.snoRNA.0004-1;Name=exon;description=exon;size=78;Parent=apidb|rna_tryp_XI-1036e06.p1k.snoRNA.0004-1
    (This is line 46748)##FASTA
    >apidb|TB927.5.300b
    ATGGCTCACGGCTCGATTCCAGTTATTGATGTCGGCCCTCTGTTCTGTGATGGAGAAAAG
    GGGATGATGGATGTTGCGAAACAGATTGATCATGCCTGTAGGACGTGGGGTGTTTTTCTT
    GTTGTGGGTCATCCCATTCCCCGTGAGCGAACGGAAAAGTTGATGGAAATGGCCAAGGCT
    TTTTTTTCGCTTCCATTGGAAGAGAAACTTAAGGTTGATATTCGAAAGAGCAAACATCAT
    CGCGGTTACGGATGCCTCGATGCGGAGAATGTTGACCCAACGAAACCATTTGATTGTAAA
    GAAACATTTAATATGGGCTGTCATCTCCCTGAGGATCACCCCGATGTTGCAGCTGGAAAG
    CCATTGCGTGGACCGAACAATCACCCCACGCAAGTGAAAGGTTGGGTAGAGTTGATGAAC
    AGACATTATCGCGAAATGCAGGAATTTGCCCTCGTTATTCTTCGTGCCCTCGCACTCGCT
    ATTGGTTTAAAGAAAGACTTTTTCGATACCAAATTTGATGAACCTTTGAGTGTGTTCCGT
    ATGCTACATTATCCTCCACAAAAGCAAGGGACCCGTTATCCCATCGTGTGTGGTGAGCAT
    ACGGATTATGGTATTATTACATTACTCTACCAAGATTCGGTGGGAGGACTGCAGGTGCGC
    AATCTGTCAGATGAGTGGGTGGATGTGGAACCCCTTGAAGGAAGTTTTGTTGTGAATATT
    GGGGACATGATGAATATGTGGAGTAATGGCCGTTACCGCTCAACACCGCATCGCGTTCGC
    TTAACCACAACTGATCGCTACTCCATGCCATTTTTCTGTCAGCCTAATCCTTATACTGTT
    ATTAAATGCCTTGATCATTGCCATTCGCCAAGCAATCCCCCCAAATATCCACCAGTCCGT
    GCTGTGGATTGGTTGCTGAAGCGTTTCGCGGAAACATATGCCCATCGCAAAACAAAGATG
    TGA
    >apidb|cds_TB927.5.300b-1
    MAHGSIPVIDVGPLFCDGEKGMMDVAKQIDHACRTWGVFLVVGHPIPRERTEKLMEMAKA
    FFSLPLEEKLKVDIRKSKHHRGYGCLDAENVDPTKPFDCKETFNMGCHLPEDHPDVAAGK
    PLRGPNNHPTQVKGWVELMNRHYREMQEFALVILRALALAIGLKKDFFDTKFDEPLSVFR
    MLHYPPQKQGTRYPIVCGEHTDYGIITLLYQDSVGGLQVRNLSDEWVDVEPLEGSFVVNI
    GDMMNMWSNGRYRSTPHRVRLTTTDRYSMPFFCQPNPYTVIKCLDHCHSPSNPPKYPPVR
    AVDWLLKRFAETYAHRKTKM

    Why does htseq-count gives now an error?
    Again, thanks for all!
    Sandra

    Comment


    • #17
      In my understanding of the GFF standard, a GFF file is supposed to contain GFF lines and not FASTA lines. I've seen this before that a full FASTA file is concatenated to the end of the GFF file but this is really confusing for any software trying to parse the file. Please remove everything that is not GFF from your GFF file.

      Comment


      • #18
        Hi SImon,

        I am really desiring now to be working with Drosophila or something similar...
        I applied what you told me:

        samtools view s_1_sequence_clipped_tophat.bam | htseq-count -s no -i ID - TbruceiTreu927_TriTrypDB-3.3.4.gff > count.txt


        Warning: Skipping read 'HWUSI-EAS582_211:1:106:1372:1342', because chromosome 'GeneDB|Tb927_01_v4', to which it has been aligned, did not appear in the GFF file

        And after I stopped the execution of the program:

        ^CError: 'itertools.chain' object has no attribute 'get_line_number_string'
        [Exception type: AttributeError, raised in count.py:198]

        Maybe it was not a good idea to eliminate the fasta part...

        Comment


        • #19
          Well, it is a only warning that told you that a single read got skipped. I suppose you can ignore that. If you have contigs in your reference to which the GFF file does not assign any exons it is completely expected to get a few of these warnings.

          Comment


          • #20
            Well, SImon, this is what I get:


            apidb|exon_Tb09.160.0220-1 0

            And over 12.000 similar, and:

            no_feature 0
            ambiguous 0
            too_low_aQual 0
            not_aligned 0
            alignment_not_unique 11678480

            But at the txt document I couldn't see anything else about the 11678480 aligments. Did it worked?? Can I see the aligments somewhere?

            And Soooo sorry for disturbing you so much. THanks for your attention!
            Sandra

            Comment


            • #21
              I would say, no, it didn't work, if you see only zeroes.

              I've just had another look at the excerpt from your GFF file that you posted above, and it does not seem right. htseq-count counts on the level of genes, not exons. Hence, for each "exon" line, the attribute "ID" (or whatever you have specified with "-i") has to be theID of the gene. All exons from the same gene must have the same ID. In your file, it seems that some exon or transcript number is appended to the ID. You may need some perl (or python or sed) script to remove these.

              Comment


              • #22
                OK, Simon...
                I will ty it...
                Thanks for your help!
                SAndra

                Comment


                • #23
                  HTSeq_ Error occured when reading first line of sam file.

                  Hello,

                  I am very new in bioinformatics and I am having some problems that can´t solve. I have seen similar threads to this one, but not exactly the same, and I am not sure how to solve my issue. When running HTSeq, I received this error:

                  100000 GFF lines processed.
                  200000 GFF lines processed.
                  217821 GFF lines processed.
                  Error occured when reading first line of sam file.

                  [Exception type: StopIteration, raised in count.py:79]
                  100000 GFF lines processed.
                  200000 GFF lines processed.
                  217821 GFF lines processed.
                  Error occured when reading first line of sam file.

                  [Exception type: StopIteration, raised in count.py:79]
                  100000 GFF lines processed.
                  200000 GFF lines processed.
                  217821 GFF lines processed.
                  Error occured when reading first line of sam file.

                  It does not tell me any specific error for why it can´t read the first line of the sam file, like in other threads that I have seen.

                  I have used Bowtie for the alignment to a reference metagenome, and samtools to create the nameSorted.sam files.

                  Any thoughts?

                  Thank you!

                  Comment


                  • #24
                    Please post the command you used and the first few lines of the SAM/BAM file in question.

                    Comment


                    • #25
                      Htseq-count throws error

                      Dear All,


                      I am using htseq count to get exonic read count fro targeted seq data..
                      I have output BAM file with alignments and I need to get the read counts for the respective exonics..the exon gtf file looks like this,
                      chr9 bed2gff exon 133589697 133589852 0 + . ABL1
                      chr9 bed2gff exon 133709410 133709441 0 + . ABL1
                      chr9 bed2gff exon 133710251 133710279 0 + . ABL1

                      And the error I receive from HT SEQ-count is,
                      Error occurred when processing GFF file (line 1 of file Regions_new.GTF):
                      Failure parsing GFF attribute line
                      [Exception type: Value Error, raised in init.py:164]
                      Pleas elet me know why I am having this error, and the command line I used is,

                      htseq-count -a 10 -m intersection-strict -s yes mydata.sam Regions_new.GTF
                      Kind regards,
                      Alva

                      Comment


                      • #26
                        That's neither GFF nor GTF format. If bed2gff produced that then the program is broken.

                        Comment


                        • #27
                          Dear Devon,
                          Thanks for the prompt reply. I would like at ask..So I have the targeted region in bed format and converted to this above mentioned format..

                          So in order to extract the read count for exons what would be most suggested tool and the format of the bed file. ??
                          Kind regards,
                          Alva

                          Comment


                          • #28
                            You're probably better off with featureCounts, since it allows counts for alignments to multiple features (this will treat spliced alignments properly).

                            Anyway, the following bit of awk code should work with your BED file:
                            Code:
                            awk '{printf("%s\tawk\texon\t%i\t%i\t.\t%s\t.\texon_id \"%s%i\"\n", $1, $2+1, $3, $6, $4, NR)}' foo.bed > foo.gff
                            Then the exons will be labeled with the associated gene and then a number (the line number). I haven't tested this, but something like that should work.

                            Comment


                            • #29
                              Ok, Thanks
                              Kind regards,
                              Alva

                              Comment


                              • #30
                                Hi everyone,

                                I have a similar problem with my gff file. I get the error:

                                Error occured when processing GFF file (line 1 of file /gpfs/scratch/cbh12wsu/Thaps3_chromosomes_geneModels_FilteredModels2.gff):
                                Feature fgenesh1_pg.C_chr_1000001 does not contain a 'gene_id' attribute
                                [Exception type: ValueError, raised in count.py:53]

                                My gff files looks like this:
                                chr_1 JGI exon 300 1153 . - . name "fgenesh1_pg.C_chr_1000001"; transcriptId 867
                                chr_1 JGI CDS 300 1153 . - 0 name "fgenesh1_pg.C_chr_1000001"; proteinId 867; exonNumber 10
                                chr_1 JGI exon 1199 2425 . - . name "fgenesh1_pg.C_chr_1000001"; transcriptId 867
                                chr_1 JGI CDS 1199 2425 . - 2 name "fgenesh1_pg.C_chr_1000001"; proteinId 867; exonNumber 9
                                chr_1 JGI exon 2512 2935 . - . name "fgenesh1_pg.C_chr_1000001"; transcriptId 867
                                chr_1 JGI CDS 2512 2935 . - 2 name "fgenesh1_pg.C_chr_1000001"; proteinId 867; exonNumber

                                I'm guessing there is something missing? Can anyone help?

                                Thanks!
                                Padmoo

                                Comment

                                Working...
                                X