[Simon Anders]

Please, read the error message. It says clearly:

"Error in newCountDataSet(countsTable, conds) :The countData is not integer."

Sorry but you are about the thousandth person to ask this, and I'm a bit desperate because I thought I wrote it nice and clearly in the manual:

DESeq expects as input a table, which tells, for each sample and each gene, how many reads fall onto the gene in the sample.

It seems that many users feel an irresistible urge to do strange things to their counts and then supply something like the number of read per million reads, the number per kilobase transcript length, or whatever.

So, please, people, stop normalizing, dividing, rounding or doing whatever with your data, and just hand over the raw unnormalized integer counts.

[sunkorner]

I am using DESeq and I get the following error.

> exampleFile = "/home/sun/Desktop/xylella_normalized_rpm.tab"
> countsTable <- read.delim(exampleFile, header=TRUE, stringsAsFactors=TRUE)
> head(countsTable)
Locus_tag WT.1 WT.2 WT.3 mutant.1 mutant.2 mutant.3
1 PD0001 32.59 35.18 200.56 37.90 173.21 224.03
2 PD0002 26.55 11.80 149.89 18.16 130.71 162.92
3 PD0003 9.73 7.46 171.14 14.50 129.91 203.61
4 PD0004 7.36 6.46 20.83 12.31 16.76 29.63
5 PD0005 47.71 69.57 481.91 82.01 372.57 567.03
6 PD0006 5.13 4.90 175.65 13.04 154.31 202.34
> rownames( countsTable) <- countsTable$Locus_tag
> countsTable <- countsTable[ ,-1 ]
> conds <- c( "W", "W", "W", "M", "M", "M" )
> cds <- newCountDataSet( countsTable, conds )
Error in newCountDataSet(countsTable, conds) :
The countData is not integer.
> cds <- newCountDataSet( countsTable, conds)
Error in newCountDataSet(countsTable, conds) :
The countData is not integer.
> head(countsTable)
WT.1 WT.2 WT.3 mutant.1 mutant.2 mutant.3
PD0001 32.59 35.18 200.56 37.90 173.21 224.03
PD0002 26.55 11.80 149.89 18.16 130.71 162.92
PD0003 9.73 7.46 171.14 14.50 129.91 203.61
PD0004 7.36 6.46 20.83 12.31 16.76 29.63
PD0005 47.71 69.57 481.91 82.01 372.57 567.03
PD0006 5.13 4.90 175.65 13.04 154.31 202.34
> str(countsTable)
'data.frame': 2123 obs. of 6 variables:
$ WT.1 : num 32.59 26.55 9.73 7.36 47.71 ...
$ WT.2 : num 35.18 11.8 7.46 6.46 69.57 ...
$ WT.3 : num 200.6 149.9 171.1 20.8 481.9 ...
$ mutant.1: num 37.9 18.2 14.5 12.3 82 ...
$ mutant.2: num 173.2 130.7 129.9 16.8 372.6 ...
$ mutant.3: num 224 162.9 203.6 29.6 567 ...

I have to still visualize my out SAM file, so I will get back on that soon.

Thank you

[Simon Anders]

htseq-count -m intersection-nonempty -o out sorted.sam chr.gff > locus_out

This won't work with your GFF file, because it does not have 'gene_id' attributes, and you don't have a '-i'.


This is why you needed this here:

htseq-count -m -t CDS -i locus_tag -o out sorted.sam chr.gff > PD_out

By the by looking into genome browser, the number of reads mapping to the gene too are counted less for my dataset. It is because some genes are overlapping each other that the parameters I gave doesnt let it to recognise. I need a count of uniquley mapped genes, partially mapped genes and not unique genes.

If you had many read falling onto overlapping genes, they would be counted as ambiguous. "no_feature" means that they don't fall anywhere.

What you should do is visualize the 'out' SAM file along your GFF file and see if you find any read that is marked as "no_feature" and appears to be on a feature. Then, make a screenshot and post it here, and we can see what is wrong.

[sunkorner]

I tried to give command,

htseq-count -m intersection-nonempty -o out sorted.sam chr.gff > locus_out

I get error line 72 gene feature doesnt exist for tRNA .......then the progam exits. When I checked the line I see that the feature in this line has 'exon'. Such lines exists for many tRNAs in this gff file. but which are other transcripts we have different information.

NC_004556.1 RefSeq gene 17033 17950 . + . ID=NC_004556.1:hemF;locus_tag=PD0015;db_xref=GeneID:1144215
NC_004556.1 RefSeq CDS 17033 17947 . + 0 ID=NC_004556.1:hemF:unknown_transcript_1;Parent=NC_004556.1:hemF;locus_tag=PD0015;EC_number=1.3.3.3;note=catalyzes the conversion of the propionic acid groups of rings I and III to vinyl groups during heme synthesis;transl_table=11;product=coproporphyrinogen III oxidase;protein_id=NP_778274.1;db_xref=GI:28197960;db_xref=GeneID:1144215;exon_number=1
NC_004556.1 RefSeq start_codon 17033 17035 . + 0 ID=NC_004556.1:hemF:unknown_transcript_1;Parent=NC_004556.1:hemF;locus_tag=PD0015;EC_number=1.3.3.3;note=catalyzes the conversion of the propionic acid groups of rings I and III to vinyl groups during heme synthesis;transl_table=11;product=coproporphyrinogen III oxidase;protein_id=NP_778274.1;db_xref=GI:28197960;db_xref=GeneID:1144215;exon_number=1
NC_004556.1 RefSeq stop_codon 17948 17950 . + 0 ID=NC_004556.1:hemF:unknown_transcript_1;Parent=NC_004556.1:hemF;locus_tag=PD0015;EC_number=1.3.3.3;note=catalyzes the conversion of the propionic acid groups of rings I and III to vinyl groups during heme synthesis;transl_table=11;product=coproporphyrinogen III oxidase;protein_id=NP_778274.1;db_xref=GI:28197960;db_xref=GeneID:1144215;exon_number=1
NC_004556.1 RefSeq gene 20728 20804 . - . locus_tag=PD0016;db_xref=GeneID:1144216
NC_004556.1 RefSeq exon 20728 20804 . - . gbkey=tRNA;locus_tag=PD0016;product=tRNA-Met;note=found by tRNAscan;db_xref=GeneID:1144216;exon_number=1

So I am exculsively giving the command htseq-count -m -t CDS -i locus_tag -o out sorted.sam chr.gff > PD_out

Then in PD_out we see,
.........
,..........
PD2111 4
PD2112 2
PD2113 7
PD2105 15
PD2106 11
PD2107 35
PD2108 10
PD2109 21
PD2110 17
.....
...
....
no_feature 14905057
ambiguous 13
too_low_aQual 0
not_aligned 1913205
alignment_not_unique 0

Something like this. Am I on the right path.

By the by looking into genome browser, the number of reads mapping to the gene too are counted less for my dataset. It is because some genes are overlapping each other that the parameters I gave doesnt let it to recognise. I need a count of uniquley mapped genes, partially mapped genes and not unique genes.

Any suggestions please.

Thank you

[sunkorner]

There are also many regions in between the genes where the reads are mapping....but the ones mapping onto the gene too are not shown in the count

[sunkorner]

Results from genome browser view,

I could see many reads actually mapping to each gene which we call locus tag. I get lesser than the actual count for each gene. Any way to increase it. When I remove the -i option then I get error, tRNAs have no gene id feature, error on 72 line and the program exit. Please let me know how to proceed.

Any input from guys working with bacterial genomes will be very helpful

Thank you

[sunkorner]

Locus tags are systematic, unique identifiers that are assigned to each gene in GenBank. Therefore I used them as identifiers in this context too. However I am not sure if the use of CDS is causing this problem of lack of identification.

I will also look into the genome browser and get back to you with results. Thank you

[Simon Anders]

Then, you should really check things with a genome browser. Check especially whether the reads are really on the same strand as the gene, to make sure that '-s yes' is correct.

[Simon Anders]

Again: Are the five lines you quote the whole output that you get from htseq-count? And have you looked at your data and the GTF file with a genome browser?

And what is a "locus tag"? (Sorry, I'm not very familiar with prokaryotic data.)

[sunkorner]

On the other hand: Is this your complete output? Not a single feature? Or did you cut them off.

I just gave the last few lines where no features is shown, I also get the count of reads mapped to each locus_tag. eg
PD0003 23
PD0001 5
etc

Thank you very much

[sunkorner]

Please find the the gff file I have used at this url ftp://ftp.ncbi.nlm.nih.gov/genomes/B.../NC_004556.gff

otherwise also please check a sample of it, hope it will be readable in the post

##gff-version 3
#!gff-spec-version 1.14
#!source-version NCBI C++ formatter 0.2
##Type DNA NC_004556.1
NC_004556.1 RefSeq source 1 2519802 . + . organism=Xylella fastidiosa Temecula1;mol_type=genomic DNA;strain=Temecula1;db_xref=taxon:183190;note=Pierce%27s disease strain
NC_004556.1 RefSeq gene 146 1465 . + . ID=NC_004556.1:dnaA;locus_tag=PD0001;db_xref=GeneID:1144259
NC_004556.1 RefSeq CDS 146 1462 . + 0 ID=NC_004556.1:dnaA:unknown_transcript_1;Parent=NC_004556.1:dnaA;locus_tag=PD0001;note=binds to the dnaA-box as an ATP-bound complex at the origin of replication during the initiation of chromosomal replication%3B can also affect transcription of multiple genes including itself.;transl_table=11;product=chromosomal replication initiation protein;protein_id=NP_778260.1;db_xref=GI:28197946;db_xref=GeneID:1144259;exon_number=1
NC_004556.1 RefSeq start_codon 146 148 . + 0 ID=NC_004556.1:dnaA:unknown_transcript_1;Parent=NC_004556.1:dnaA;locus_tag=PD0001;note=binds to the dnaA-box as an ATP-bound complex at the origin of replication during the initiation of chromosomal replication%3B can also affect transcription of multiple genes including itself.;transl_table=11;product=chromosomal replication initiation protein;protein_id=NP_778260.1;db_xref=GI:28197946;db_xref=GeneID:1144259;exon_number=1
NC_004556.1 RefSeq stop_codon 1463 1465 . + 0 ID=NC_004556.1:dnaA:unknown_transcript_1;Parent=NC_004556.1:dnaA;locus_tag=PD0001;note=binds to the dnaA-box as an ATP-bound complex at the origin of replication during the initiation of chromosomal replication%3B can also affect transcription of multiple genes including itself.;transl_table=11;product=chromosomal replication initiation protein;protein_id=NP_778260.1;db_xref=GI:28197946;db_xref=GeneID:1144259;exon_number=1
NC_004556.1 RefSeq gene 1747 2847 . + . ID=NC_004556.1:dnaN;locus_tag=PD0002;db_xref=GeneID:1144260
NC_004556.1 RefSeq CDS 1747 2844 . + 0 ID=NC_004556.1:dnaN:unknown_transcript_1;Parent=NC_004556.1:dnaN;locus_tag=PD0002;EC_number=2.7.7.7;note=binds the polymerase to DNA and acts as a sliding clamp;transl_table=11;product=DNA polymerase III subunit beta;protein_id=NP_778261.1;db_xref=GI:28197947;db_xref=GeneID:1144260;exon_number=1
NC_004556.1 RefSeq start_codon 1747 1749 . + 0 ID=NC_004556.1:dnaN:unknown_transcript_1;Parent=NC_004556.1:dnaN;locus_tag=PD0002;EC_number=2.7.7.7;note=binds the polymerase to DNA and acts as a sliding clamp;transl_table=11;product=DNA polymerase III subunit beta;protein_id=NP_778261.1;db_xref=GI:28197947;db_xref=GeneID:1144260;exon_number=1
NC_004556.1 RefSeq stop_codon 2845 2847 . + 0 ID=NC_004556.1:dnaN:unknown_transcript_1;Parent=NC_004556.1:dnaN;locus_tag=PD0002;EC_number=2.7.7.7;note=binds the polymerase to DNA and acts as a sliding clamp;transl_table=11;product=DNA polymerase III subunit beta;protein_id=NP_778261.1;db_xref=GI:28197947;db_xref=GeneID:1144260;exon_number=1
NC_004556.1 RefSeq gene 3153 4247 . + . ID=NC_004556.1:recF;locus_tag=PD0003;gene_synonym=uvrF;db_xref=GeneID:1144261
NC_004556.1 RefSeq CDS 3153 4244 . + 0 ID=NC_004556.1:recF:unknown_transcript_1;Parent=NC_004556.1:recF;locus_tag=PD0003;gene_synonym=uvrF;note=Required for DNA replication%3B binds preferentially to single-stranded%2C linear DNA;transl_table=11;product=recombination protein F;protein_id=NP_778262.1;db_xref=GI:28197948;db_xref=GeneID:1144261;exon_number=1
NC_004556.1 RefSeq start_codon 3153 3155 . + 0 ID=NC_004556.1:recF:unknown_transcript_1;Parent=NC_004556.1:recF;locus_tag=PD0003;gene_synonym=uvrF;note=Required for DNA replication%3B binds preferentially to single-stranded%2C linear DNA;transl_table=11;product=recombination protein F;protein_id=NP_778262.1;db_xref=GI:28197948;db_xref=GeneID:1144261;exon_number=1
NC_004556.1 RefSeq stop_codon 4245 4247 . + 0 ID=NC_004556.1:recF:unknown_transcript_1;Parent=NC_004556.1:recF;locus_tag=PD0003;gene_synonym=uvrF;note=Required for DNA replication%3B binds preferentially to single-stranded%2C linear DNA;transl_table=11;product=recombination protein F;protein_id=NP_778262.1;db_xref=GI:28197948;db_xref=GeneID:1144261;exon_number=1

I am looking for the count of all the reads mapped to each gene. Please let me know if it is correct way I am doing using locus_tag and CDs. It is also microbial genome.

Thank you

[Simon Anders]

By the by what is the no_feature mean: are these number of reads not mapped?

No, it means that the read does not overlap with any feature in your GTF file, i.e., that it falls onto intergenic or intronic space.

To investigate, first display your SAM file and you GTF file alongside in a genome browser (e.g., IGV) and inspect visually, whether there are indeed many reads that fall between features.

On the other hand: Is this your complete output? Not a single feature? Or did you cut them off.

If not, are you sure your GFF file fits your options? You have instructed htseq-count to only use those lines in your GFF file that have the string "CDS" in its third field, and have a suitable attribute called "locus_tag" in the last field. Maybe post an extract of the GFF file.

[sunkorner]

mbobro2, may be do u think we need to sort the sam file before HTSeq count, sometimes if we dont sort it can give odd results

[sunkorner]

I am using the following commands, I get lot of reads as "no_feature", I tried union as well as intersection strict, but no much change in the nubmer of no_feature. only the ambiguous reads nubmer has change. My total number of reads is around 18 million. So I am worried in this regard.

By the by what is the no_feature mean: are these number of reads not mapped?

Please suggest what might be wrong with my htseq command. Thank you

===
htseq-count -m union -s yes -t CDS -i locus_tag -o WT-1_geneid_htseq_out3 sorted.sam /home/sun/Rna_seq_data/Xylella/chrComplete.gff > PD_out_union &
no_feature 14867530
ambiguous 4453
too_low_aQual 0
not_aligned 1913205
alignment_not_unique 0

htseq-count -m intersection-strict -s yes -t CDS -i locus_tag -o WT-1_geneid_htseq_out2 sorted.sam /home/sun/Rna_seq_data/Xylella/chrComplete.gff > PD_out &
no_feature 14905057
ambiguous 13
too_low_aQual 0
not_aligned 1913205
alignment_not_unique 0

[Simon Anders]

Have you looked into you SAM file at the locus of the gene with a gene browser (e.g., IGV)?