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 12-12-2011, 04:10 AM   #21
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

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.
Simon Anders is offline   Reply With Quote
Old 12-12-2011, 06:25 AM   #22
Dedeusan
Junior Member
 
Location: spain

Join Date: Nov 2011
Posts: 7
Default

OK, Simon...
I will ty it...
Thanks for your help!
SAndra
Dedeusan is offline   Reply With Quote
Old 01-21-2014, 07:36 AM   #23
m232
Junior Member
 
Location: woods hole

Join Date: Aug 2013
Posts: 2
Unhappy 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!
m232 is offline   Reply With Quote
Old 01-21-2014, 07:40 AM   #24
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Please post the command you used and the first few lines of the SAM/BAM file in question.
dpryan is offline   Reply With Quote
Old 03-26-2015, 07:11 AM   #25
alvarani
Junior Member
 
Location: Sweden

Join Date: Mar 2013
Posts: 5
Default 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
alvarani is offline   Reply With Quote
Old 03-26-2015, 07:14 AM   #26
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

That's neither GFF nor GTF format. If bed2gff produced that then the program is broken.
dpryan is offline   Reply With Quote
Old 03-26-2015, 07:23 AM   #27
alvarani
Junior Member
 
Location: Sweden

Join Date: Mar 2013
Posts: 5
Default

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
alvarani is offline   Reply With Quote
Old 03-26-2015, 07:35 AM   #28
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 03-27-2015, 02:57 AM   #29
alvarani
Junior Member
 
Location: Sweden

Join Date: Mar 2013
Posts: 5
Smile

Ok, Thanks
__________________
Kind regards,
Alva
alvarani is offline   Reply With Quote
Old 06-19-2015, 09:14 AM   #30
padmoo
Member
 
Location: Norwich

Join Date: Jun 2015
Posts: 16
Default

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
padmoo is offline   Reply With Quote
Old 06-19-2015, 09:40 AM   #31
cmbetts
Senior Member
 
Location: Bay Area

Join Date: Jun 2012
Posts: 118
Default

Quote:
Originally Posted by padmoo View Post
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
In the gff format, the final column has a semi colon seperated list of attributes describing the feature. htseq-counts is looking for one called gene_id that lets it assign counts on the gene level, but your gff doesn't have that as one of the attributes. I'm not sure if "fgenesh1_pg.C_chr_1000001" is supposed to be a gene of if it's something else (it would make your life much easier if it is), but somehow you need to add something that looks like 'gene_id "myGene"' to the gff where it describes exons.
cmbetts is offline   Reply With Quote
Old 06-19-2015, 10:01 AM   #32
padmoo
Member
 
Location: Norwich

Join Date: Jun 2015
Posts: 16
Default

Hi cmbetts,

thanks for replying so quickly. I tried to find some other gff files and came across one that looks like this:
##gff-version 3
##sequence-region bd_19x4 1 66660
bd_19x4 prediction gene 639 1694 0 - . ID=gene_1;Name=jgi.p|Thaps3_bd|833;portal_id=Thaps3_bd;proteinId=833;transcriptId=833
bd_19x4 prediction mRNA 639 1694 . - . ID=mRNA_1;Name=jgi.p|Thaps3_bd|833;Parent=gene_1;proteinId=833;track=FilteredModels1;transcriptId=833
bd_19x4 prediction exon 639 1694 . - . ID=exon_1_1;Parent=mRNA_1
bd_19x4 prediction CDS 639 1694 . - 0 ID=CDS_1;Parent=mRNA_1
bd_19x4 prediction gene 3895 5607 0 - . ID=gene_2;Name=jgi.p|Thaps3_bd|1845;portal_id=Thaps3_bd;proteinId=1845;transcriptId=1845
bd_19x4 prediction mRNA 3895 5607 . - . ID=mRNA_2;Name=jgi.p|Thaps3_bd|1845;Parent=gene_2;proteinId=1845;track=FilteredModels1;transcriptId=1845
bd_19x4 prediction exon 5551 5607 . - . ID=exon_2_1;Parent=mRNA_2
bd_19x4 prediction five_prime_UTR 5551 5607 . - . ID=UTR5_2;Parent=mRNA_2
bd_19x4 prediction exon 3895 5455 . - . ID=exon_2_2;Parent=mRNA_2
bd_19x4 prediction five_prime_UTR 5097 5455 . - . ID=UTR5_2;Parent=mRNA_2

This one has an ID=gene but it's only numbered. Maybe there are no gene models?
I'm really confused... The files are from the JGI genome portal.
padmoo is offline   Reply With Quote
Old 06-20-2015, 07:36 AM   #33
padmoo
Member
 
Location: Norwich

Join Date: Jun 2015
Posts: 16
Default

Solved my problem. Didn't pay attention to the -i option...
padmoo is offline   Reply With Quote
Old 06-21-2015, 07:39 AM   #34
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Be careful here. In this file, the "exon" lines only contain transcript IDs ("mRNA_1") but not gene IDs ("gene_1") in their "Parent" fields. This will not work, because if an exon appears in several transcripts of the same gene, it will be listed in several lines in the GFF file, each with a different "Parent" value, and if you use "-i Parent", htseq-count will think that these different exon lines correspond to different genes and consider the assignment "ambiguous" rather than counting it for the gene.

You will have to tweak the GFF3 file such that you add to each "exon" line an additional attribute, say, "gene", which contains the gene ID rather than the transcript of mRNA ID, i.e. (in the logic of GFF3), the Parent ID of the Parent.

Sorry that this is so complicated, but unfortunately, these file formats are all ill-specified and keep changing. (Maybe I am getting senile, but I am pretty sure that five years ago, when I wrote htseq-count, the GFF3 specification (or was it maybe only GFF2 then?) had "exon" and "gene" lines but no "mRNA" lines, and the "Parent" of an "exon" line was a "gene" line.
Simon Anders is offline   Reply With Quote
Old 07-20-2015, 07:09 AM   #35
populus
Junior Member
 
Location: Beijing

Join Date: Jul 2015
Posts: 3
Default

Hi simon,

Recently, I am trying to do DEG analysis. I use HTSeq-count to count reads. Here I want to know do I need use the .gtf file during reads mapping with Tophat (-G) before using the HTSeq-count? Because I found there was a little difference between the output.counts. Thank you very much.
populus is offline   Reply With Quote
Old 07-20-2015, 10:38 PM   #36
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

While you don't need to supply the GTF file to tophat, it will often improve the results slightly.
dpryan is offline   Reply With Quote
Old 07-21-2015, 07:53 AM   #37
populus
Junior Member
 
Location: Beijing

Join Date: Jul 2015
Posts: 3
Default

Quote:
Originally Posted by dpryan View Post
While you don't need to supply the GTF file to tophat, it will often improve the results slightly.
Thank you for the reply. I have another question. When I used HTSeq-count to count the read pairs for paired-end data, why I just get the counts of reads but not the read pairs. e.g. for my .sam, it has 15241301 lines excluded headers. After HTSeq-count(-r name -s no) analysis, I got the output.counts:
(sum of counts of genes: 14939019)
__no_feature 129677
__ambiguous 172605
__too_low_aQual 0
__not_aligned 0
__alignment_not_unique 0

So how can I get the read pairs? Thank you very much.

ps: After tophat, I samtools rmdup the accepted_hits.bam, filtered the unique_mapping_reads and sorted the .bam.

Last edited by populus; 07-21-2015 at 08:09 AM. Reason: add information
populus is offline   Reply With Quote
Old 07-21-2015, 10:54 PM   #38
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

1) Don't remove duplicates from RNAseq data without a very very good reason.
2) HTseq-count is giving you counts of read pairs, or at least it should be if you actually aligned things as pairs.
dpryan is offline   Reply With Quote
Old 07-22-2015, 06:00 AM   #39
populus
Junior Member
 
Location: Beijing

Join Date: Jul 2015
Posts: 3
Default

Quote:
Originally Posted by dpryan View Post
1) Don't remove duplicates from RNAseq data without a very very good reason.
2) HTseq-count is giving you counts of read pairs, or at least it should be if you actually aligned things as pairs.
I think I have figured out the problem. I found there were suffix .1 and .2 in the seq_name of paired reads. It must be counted as single end data. I removed the suffix and the htseq-count worked well.

Thank you very much.
populus is offline   Reply With Quote
Old 07-20-2016, 06:03 PM   #40
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default Most of my reads are going to no_feature

Good evening SEQanswers community,

I am having issues with counting my reads using HTseq. I found this thread, but I was not able to use the same modifications to help my situation. First I sorted my unsorted mapped reads sam file (not sure if this was necessary), then I edited the sorted file, then I edited the gff3 file, and executed the HTSeq.scripts.count. Most of my reads went to no_feature. What am I doing wrong? Please see my input and sections of my output.

I would greatly appreciate your help!

#Sorting my mapped reads with:
samtools sort -n -O sam -T Project.sort -o Project.sort.sam Project_mapped.sam

Output:
@SQ SN:lcl|NC_017304.1_cds_WP_003516470.1_3012 LN:1380
@SQ SN:lcl|NC_017304.1_cds_WP_003516468.1_3013 LN:621
@SQ SN:lcl|NC_017304.1_cds_WP_003513356.1_3014 LN:876
@SQ SN:lcl|NC_017304.1_cds_WP_003513354.1_3015 LN:216
@SQ SN:lcl|NC_017304.1_cds_WP_003513352.1_3016 LN:381
@SQ SN:lcl|NC_017304.1_cds_WP_003513349.1_3017 LN:135
@PG ID:bowtie2 PN:bowtie2 VN:2.2.9 CL:"<removed> basic-0 --threads 12 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x Project_b2index -S Project
_mapped.sam -U /tmp/15217.unp"
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1236:29817 4 * 0 0 * * 0 0 AACAACAAAAGACAAAAAACTTAAAAACAAAAGACAAAGTCCAAAATTA BBBB/FFFFFFBBB
F/FFFB/BFFFFF/FFBFF<<<<<F/FBFFF<F<F YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1301:27869 4 * 0 0 * * 0 0 TAAGAGTTTTTCGAGGGAATAAATAACAACGTTCAAGGAAATCCT BBB<<B//<B//<<////</<B/FFFB////<FB<//</<BB/FB YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1302:11737 0 lcl|NC_017304.1_cds_WP_003514647.1_448 142 42 48M * 0 0 AAAGCCATAGATGCAGCGGTTAATGACATGGCTCTGATAAGTGGTCAG BBBBBFFBF////BFFF</F/FFF/FFFFFFFB/<<FFBFBF<FFFFB AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:11A36 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1318:93081 0 lcl|NC_017304.1_cds_WP_003514647.1_448 321 40 44M * 0 0 CCTGCCAAGAGTAAGTGACTTTTGAGGTGTTCCTGTAAATTCTT BBBBBF///</<B<///<</<</<<<//<<BFFF/<</FBBFFB AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:15G6A21 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1324:19877 0 lcl|NC_017304.1_cds_WP_014522594.1_623 2717 1 50M * 0 0 GCGGTATACCATCCAAGGGAATAGCAAACTGTGACTTTGTATACAGCTAT BBBBBBFFFFFFF/FB<//BFFFFFFFFFBFF<FF/FFFFFFFFFFBFFF AS:i:0 XS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:50 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1325:83438 0 lcl|NC_017304.1_cds_WP_003517434.1_1009 2051 24 41M * 0 0 TCAACACTTCAATCAATGCGGTTATTGAGCTGGTTAATTCC B<BBBBBF/</<B<F/<//BB//7F/<F<F/</<FFFF<FB AS:i:-9 XN:i:0 XM:i:3 XO:i:0 XG:i:0 NM:i:3 MD:Z:8G6G9G15 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1331:65339 0 lcl|NC_017304.1_cds_WP_003513122.1_105 10 40 43M * 0 0 GTTACTTTAAACCTAAGAGTTGAGGCAATAAAAAACAACGAAG <BBBFF/<//<F////</<<//<FBF7FFF//<FFB//</<FF AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:6A6A29 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1345:80429 0 lcl|NC_017304.1_cds_WP_003512596.1_1851 207 42 49M * 0 0 TGTTGAAGAGACGGGACTTCCTATCTGCCTGCATCTTG

#Edited the sorted file with sed
cp ./Project.sort.sam ./Project.sort.sed.sam
sed -i 's/lcl|NC............................_/gene/g' Project.sort.sed.sam

Pre-edit:
@SQ SN:lcl|NC_017304.1_cds_WP_003516470.1_3012 LN:1380
@SQ SN:lcl|NC_017304.1_cds_WP_003516468.1_3013 LN:621
@SQ SN:lcl|NC_017304.1_cds_WP_003513356.1_3014 LN:876
@SQ SN:lcl|NC_017304.1_cds_WP_003513354.1_3015 LN:216
@SQ SN:lcl|NC_017304.1_cds_WP_003513352.1_3016 LN:381
@SQ SN:lcl|NC_017304.1_cds_WP_003513349.1_3017 LN:135
@PG ID:bowtie2 PN:bowtie2 VN:2.2.9 CL:<removed> basic-0 --threads 12 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x cipA_b2index -S cipA
_mapped_DS3_3-5.sam -U /tmp/15217.unp"
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1236:29817 4 * 0 0 * * 0 0 AACAACAAAAGACAAAAAACTTAAAAACAAAAGACAAAGTCCAAAATTA BBBB/FFFFFFBBB
F/FFFB/BFFFFF/FFBFF<<<<<F/FBFFF<F<F YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1301:27869 4 * 0 0 * * 0 0 TAAGAGTTTTTCGAGGGAATAAATAACAACGTTCAAGGAAATCCT BBB<<B//<B//<<////</<B
/FFFB////<FB<//</<BB/FB YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1302:11737 0 lcl|NC_017304.1_cds_WP_003514647.1_448 142 42 48M * 0 0 AAAGCCATAGATGCAGCGGTTAATGACATGGCTCTGAT
AAGTGGTCAG BBBBBFFBF////BFFF</F/FFF/FFFFFFFB/<<FFBFBF<FFFFB AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:11A36 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1318:93081 0 lcl|NC_017304.1_cds_WP_003514647.1_448 321 40 44M * 0 0 CCTGCCAAGAGTAAGTGACTTTTGAGGTGTTCCTGTAA
ATTCTT BBBBBF///</<B<///<</<</<<<//<<BFFF/<</FBBFFB AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:15G6A21 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1324:19877 0 lcl|NC_017304.1_cds_WP_014522594.1_623 2717 1 50M * 0 0 GCGGTATACCATCCAAGGGAATAGCAAACTGTGACTTTGTATACAGCTAT BBBBBBFFFFFFF/FB<//BFFFFFFFFFBFF<FF/FFFFFFFFFFBFFF AS:i:0 XS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:50 YT:Z:UU

Post-edit:
@SQ SN:gene3013 LN:621
@SQ SN:gene3014 LN:876
@SQ SN:gene3015 LN:216
@SQ SN:gene3016 LN:381
@SQ SN:gene3017 LN:135
@PG ID:bowtie2 PN:bowtie2 VN:2.2.9 CL:<removed> basic-0 --threads 12 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x cipA_b2index -S cipA
_mapped_DS3_3-5.sam -U /tmp/15217.unp"
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1236:29817 4 * 0 0 * * 0 0 AACAACAAAAGACAAAAAACTTAAAAACAAAAGACAAAGTCCAAAATTA BBBB/FFFFFFBBB
F/FFFB/BFFFFF/FFBFF<<<<<F/FBFFF<F<F YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1301:27869 4 * 0 0 * * 0 0 TAAGAGTTTTTCGAGGGAATAAATAACAACGTTCAAGGAAATCCT BBB<<B//<B//<<////</<B
/FFFB////<FB<//</<BB/FB YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1302:11737 0 gene448 142 42 48M * 0 0 AAAGCCATAGATGCAGCGGTTAATGACATGGCTCTGATAAGTGGTCAG BBBBBFFBF////B
FFF</F/FFF/FFFFFFFB/<<FFBFBF<FFFFB AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:11A36 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1318:93081 0 gene448 321 40 44M * 0 0 CCTGCCAAGAGTAAGTGACTTTTGAGGTGTTCCTGTAAATTCTT BBBBBF///</<B<///<</<<
/<<<//<<BFFF/<</FBBFFB AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:15G6A21 YT:Z:UU
HISEQ-WALDORF:257:C95T2ANXX:2:1101:1324:19877 0 gene623 2717 1 50M * 0 0 GCGGTATACCATCCAAGGGAATAGCAAACTGTGACTTTGTATACAGCTAT BBBBBBFFFFFFF/
FB<//BFFFFFFFFFBFF<FF/FFFFFFFFFFBFFF AS:i:0 XS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:50 YT:Z:UU


#Edited the GFF3 file with Excel, then used ren *.txt *.gff to convert back.

Pre-edit:
##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build ASM18492v1
#!genome-build-accession NCBI_Assembly:GCF_000184925.1
##sequence-region NC_017304.1 1 3561619
##species https://www.ncbi.nlm.nih.gov/Taxonom....cgi?id=637887
NC_017304.1 RefSeq region 1 3561619 . + . ID=id0;Dbxref=taxon:637887;Is_circular=true;Name=ANONYMOUS;gbkey=Src;genome=chromosome;mol_type=genomic DNA;strain=DSM 1313
NC_017304.1 RefSeq gene 212 1543 . + . ID=gene0;Name=CLO1313_RS00010;gbkey=Gene;gene_biotype=protein_coding;locus_tag=CLO1313_RS00010;old_locus_tag=Clo1313_0001
NC_017304.1 Protein Homology CDS 212 1543 . + 0 ID=cds0;Parent=gene0;Dbxref=Genbank:WP_003516465.1;Name=WP_003516465.1;gbkey=CDS;product=chromosomal replication initiator protein DnaA;protein_id=WP_003516465.1;transl_table=11
NC_017304.1 RefSeq gene 1793 2893 . + . ID=gene1;Name=CLO1313_RS00015;gbkey=Gene;gene_biotype=protein_coding;locus_tag=CLO1313_RS00015;old_locus_tag=Clo1313_0002
NC_017304.1 Protein Homology CDS 1793 2893 . + 0 ID=cds1;Parent=gene1;Dbxref=Genbank:WP_003513339.1;Name=WP_003513339.1;gbkey=CDS;product=DNA polymerase III subunit beta;protein_id=WP_003513339.1;transl_table=11
NC_017304.1 RefSeq gene 2909 3115 . + . ID=gene2;Name=CLO1313_RS00020;gbkey=Gene;gene_biotype=protein_coding;locus_tag=CLO1313_RS00020;old_locus_tag=Clo1313_0003

Post-edit:
##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build ASM18492v1
#!genome-build-accession NCBI_Assembly:GCF_000184925.1
##sequence-region NC_017304.1 1 3561619
##species https://www.ncbi.nlm.nih.gov/Taxonom....cgi?id=637887
NC_017304.1 RefSeq region 1 3561619 . + . ID=id0
NC_017304.1 RefSeq gene 212 1543 . + . ID=gene0
NC_017304.1 Protein Homology CDS 212 1543 . + 0 ID=cds0
NC_017304.1 RefSeq gene 1793 2893 . + . ID=gene1
NC_017304.1 Protein Homology CDS 1793 2893 . + 0 ID=cds1
NC_017304.1 RefSeq gene 2909 3115 . + . ID=gene2
NC_017304.1 Protein Homology CDS 2909 3115 . + 0 ID=cds2
NC_017304.1 RefSeq gene 3133 4242 . + . ID=gene3
NC_017304.1 Protein Homology CDS 3133 4242 . + 0 ID=cds3
NC_017304.1 RefSeq gene 4410 4724 . + . ID=gene4
NC_017304.1 Protein Homology CDS 4410 4724 . + 0 ID=cds4


#Then did the count
#Count reads using htseq-count
python2.7 -m HTSeq.scripts.count -t gene -i ID Project.sort.sed.sam GCF_000184925.1_2.gff > Project_counts.txt

Output: Project_counts.txt

gene995^M 0
gene996^M 0
gene997^M 0
gene998^M 0
gene999^M 0
__no_feature 14821716
__ambiguous 0
__too_low_aQual 999116
__not_aligned 6400468
__alignment_not_unique 0

Last edited by jmwhitha; 07-20-2016 at 06:05 PM. Reason: need to know
jmwhitha 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 09:40 AM.


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