SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
DEXSeq Using Counts File From htseq-count FuzzyCoder Bioinformatics 20 01-03-2016 11:18 PM
Error with GTF file when using htseq-count MDonlin Bioinformatics 13 01-13-2015 08:29 AM
Strange error when using htseq-count shhuang Bioinformatics 13 11-19-2012 12:40 AM
htseq count missing last mapped fragment mapardo RNA Sequencing 2 11-08-2011 05:27 AM
Htseq-count Vs CountOverlap function in IRanges biofreak General 5 06-29-2011 10:36 AM

Reply
 
Thread Tools
Old 07-20-2016, 06:12 PM   #41
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

Does it have to do with my S? http://www-huber.embl.de/users/ander...doc/count.html says that "__no_feature: reads (or read pairs) which could not be assigned to any feature (set S as described above was empty)." I know there are start and stop positions for the GFF3. Are there corresponding start and stop positions in the sam file? If so, what column? I did look at http://genome.sph.umich.edu/wiki/SAM, but I still can't figure out if there is corresponding start and stop info. Thanks again for the help!

http://www-huber.embl.de/users/ander...doc/count.html
jmwhitha is offline   Reply With Quote
Old 07-20-2016, 11:57 PM   #42
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Quote:
@SQ SN:gene3013 LN:621
So, each contig/chromosome in your SAM file is a gene. This looks like you aligned to a transcriptome but then tried to proceed with a work flow meant for SAM files aligned to a genome.

I think you got confused about what a GFF3 file is.

Better explain a bit more what you are trying to do.
Simon Anders is offline   Reply With Quote
Old 07-21-2016, 06:12 AM   #43
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default Multi-fasta

Thank you very much for your fast reply! Yes, each contig/chromosome (aka @SQ) is actually a gene. The format of the multi-fasta file I used for my bowtie2 alignment is below. It's not a transcriptome (http://www.ncbi.nlm.nih.gov/nuccore/NC_017304). Why do you think that? My goal is to map count gene reads, not transcript reads. I thought by replacing "lcl|NC............................_" with gene, I could get HTseq count to count my reads. What should I do instead? Thank you again!

>lcl|NC_017304.1_cds_WP_003516465.1_1 [gene=CLO1313_RS00010] [protein=/inference=EXISTENCE: similar to AA sequence:SwissProt:A3DHZ4.1] [protein_id=WP_003516465.1] [location=212..1543
]
ATGAATACTCAGTTGAATGAAATATGGCAAAAAACTTTAGGACTGCTTAAAAATGAGCTTACAGAAATCA
GTTTTAACACCTGGATCAAGACCATCGATCCATTGTCCTTGACAGGCAATACTATAAACCTGGCTGTCCC
GGCGGAATTCAACAAGGGAATTCTTGAGTCCAGGTATCAAACTCTGATTAAAAATGCCATTAAGCAAGTT
ACTTTTAAGGAATACGAGATTGCATTTATCGTGCCTTCACAGGAAAATTTAAACAAGCTGACGAAGCAGA
CCGAGTCCGCCGGCAATGAGGATTCTCCTTTGTCAGTATTAAACCCGAAGTACACGTTTGACACTTTTGT
CATAGGAAACAGCAACAGATTTGCACACGCAGCCGCACTGGCCGTGGCCGAGGCACCGGGAAAAGCATAC
AATCCCTTGTTCATATATGGCGGAGTGGGACTTGGGAAGACTCATCTTATGCATGCCATCGGGCACTACA
TTCTGGAACAGAATTCTTCCCAAAGGGTTTTGTATGTTTCATCTGAAAAATTTACCAACGAACTTATCAA
TGCCATTAAAGACAACAGAAATGAAGAATTCAGATCCAAATACAGAAATATTGACGTACTGCTTATAGAC
GACATACAATTCATTGCCGGAAAGGAAAGAACGGAGGAGGAGTTCTTCCATACCTTCAATGCTCTTTACG
AAGCAAACAAACAGATAATCCTGTCAAGCGACAAGCCTCCGAAAGAAATTTCTCTTGAGGACCGCCTGAG
ATCCAGGTTTGAATGGGGCTTGATTGCGGACATGCAGGCACCGGATCTGGAAACCAGGATAGCAATACTA
AGGAAAAAAGCCCAGCTTGAAAACCTTACTGTTCCAAATGAAGTAATTGTATTCATTGCAGACAAGATAG
CATCAAACATCAGAGAACTTGAAGGTGCCTTAAACAGAGTAATAGCATATTCATCGCTTACGGAAAACGA
AATTACCGTCGAACTCGCCAGCGAAGCATTAAAAGACATACTGTCAGCAAACAAGGCGAAAGTTTTAAAC
TGCACCACAATCCAGGAAGCAGTGGCCAGATACTTTGACATAAGACCGGAAGAATTTAAATCAAAGAAGA
GGACAAGGGACATCGCATTCCCAAGACAAATTGCAATGTACCTGTGCAGAGAACTTACCGAAATGTCCCT
CCCAAAAATCGGCGAGGAATTCGGCGGAAGAGATCATACTACTGTAATACATGCATGTGAAAAGATAAGT
GAAGAAATCGAAAGCAACTCCGAAACCAGGAGGGCCGTAAGTGAAATAAAGAGGAACCTGCTGGGAAAAT
AA
>lcl|NC_017304.1_cds_WP_003513339.1_2 [gene=CLO1313_RS00015] [protein=/inference=EXISTENCE: similar to AA sequence:RefSeq:WP_003513339.1] [protein_id=WP_003513339.1] [location=1793..
2893]
ATGAAAATAGTTTGTTCCAAAGAACAGCTAATGGAAGGAATCAACGTCGTGCAAAAAGCAGTGCCGACAA
AAGCCACTCTAACCATACTGGAAGGAATATTGCTGGAAGCATACGACAATTTTAAAATGACCGGAAATGA
TTTGGAACTGGGAATAGAATGCCTTATAGATGCAGACATTCTGGAAAAAGGATCTATAGTCTTAAATTCA
AAAATGTTCGGAGACATAGTAAGAAGACTTCCCGACTCAGAGGTACTTATTGAAGTTAAAGAGAACAATA
CAGTTATCATTGAATGTGACAACTCTCACTTTGAGTTAAGGGGTATGCCTTCTGACAGCTTTCCGTCACT
GCCTTCAATTGAAAAAGAGAACATGATCAAAGTCAGCCAAAAGGCAATCAGGGATATGATAAGACAAACA
CTTTTTGCCGTAAGTATGGAAGGAACCAGACCGATACTTACCGGTTCACTTATTGAATGTGCAGGAAACG
AAATTACCTTCGTTTCAATAGACGGATTCAGAATGGCTCTGAGAAAAAACTTTAACAACGAAGGATTTTC
CGAATTCAGTGTTGTCGTACCCGCAAAAACCCTCAGCGAGATAGGCAAAATCTTACAGCCGGTTGATGAA
GATATTTACATATACAGTTCTCAAAACCAGATACTGTTTGAAATTGGAAATTGCAAAGTTGTATCAAGAC
TTTTAGAGGGTGAATATCTAAACTATAAAAGTATTATACCACCGGAATATGAAACCAGCGTAAGACTTAG
AACCGAGGACCTTTTGTCCAGCCTTGAAAGGGCGTCATTGATTACTTCGGACGAAAAGAAATACCCGGTT
AAATTTAATATTATAGACGATAAAATCATAATTACCTCCAACACTGAAATAGGAGCAGTAAGGGAAGAAA
TCAGAGTCGAAGTAAACGGCAGCAACATGGAAGTGGGCTTCAACCCCAGATATTTTATCGAAGCGCTCAG
GGTCATAGATGACGAGCTGGTTGACATATACTTCAATTCAAGTGTCGGTCCGTGTACAATAAGACCTCTT
GAAGGCGACAGTTTTGCATACATGATACTTCCGGTAAGAATAAATAAATAA

Commands:

#Index
bowtie2-build -o 2 --threads 12 -f NC_017304.1.nucleotide.fa cipA_b2index

#Mapped 3 fastq files to the index
bowtie2 --threads 12 --very-sensitive-local -x cipA_b2index -U 1.fastq.gz,2.fastq.gz,3.fastq.gz -S Project_mapped.sam
jmwhitha is offline   Reply With Quote
Old 07-23-2016, 09:55 AM   #44
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

If each gene is on its own pseudo-contig, you don't need htseq-count. You just count how often you see each gene ID in the third column of your SAM file which contains the chromosome to which the read was mapped. You may want to use a suitable script to skip over multi-mapping reads, though (in the easiest case, grep for all reads with NH:i:1 or something like that).

On the other hand, it is not a good idea to use a tool chain in a manner not intended by the developers of the tools unless you know enough about the tools\ internals to be sure that this is sound. Bowtie is not meant to map to a reference made up not of complete chromosomes of contigs but of individual genes, and htseq-count is neither meant to be used in this manner.

However, there are now tools designed for your manner of operation, most notably salmon and kallisto. So, the easiest and cleanest solution would be to use one of these two.

Last edited by Simon Anders; 07-23-2016 at 09:57 AM.
Simon Anders is offline   Reply With Quote
Old 08-20-2016, 09:21 AM   #45
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

Simon,

Thanks again for your quick replies. I just wanted to mention that I really appreciate your last response. While you may not be interested, I wanted to describe for others at least how I was not using the HtSeq-count tool/workflow properly. Instead of using the genome fasta file, I was using a multi-fasta file (cds only). Once I used the genome fasta file, almost everything worked properly. The one issue I ran into was that HtSeq count did not recognize the values in the first column of the GFF3 file. The reason was because the value in the first column was different from the sam file headers. The way I solved this was by changing the genome fasta file header to value in the first column of the GFF3 file before mapping.

Thank you and God bless,
Jason
jmwhitha is offline   Reply With Quote
Old 11-04-2016, 03:32 PM   #46
axa9070
Member
 
Location: Rochester, NY

Join Date: Oct 2014
Posts: 11
Default

Quote:
Originally Posted by Simon Anders View Post
Try '-i ID', not '-i ID=gene'. Only the part before the '=' is the attribute name.
For me it was Name, but this thread was very helpful.
axa9070 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:11 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO