dexseq-count.py DEXSeq gives no results

[annaprotasio]

Hello,

I am running into a bit of a problem with dexseq-count.py from DEXSeq.

dexseq-count.py. was installed and it runs OK with the sorted bam file and the gtf I provide (see below) but there are no counts in the output - not even the exons are annotated in the output. My output file looks like:

_ambiguous	0
_ambiguous_readpair_position	290728
_empty	12569972
_lowaqual	0
_notaligned	0

My sequencing data is paired and the bam file has been sorted by name (samtools sort -n file.bam). This is how it looks like:

HS21_13853:1:1101:1000:39323#6	99	Schisto_mansoni.Chr_ZW	29878411	50	97M3243N3M	=	29878440	3372	NGCCGAACCAACAGATAATTTGAACTAATAATGGATCAAATGAATAACTCAAATCAGCTGTATCCGAATAATAATACTGAATTAGTGATTTTTTATTTTT	!4=DDDFFHHHHHJIJJJJJJJHHIIIIIJJJJJFHJJJJJJIIIIJJJJJJJJJJJJJJHIJJJJHHHHHHHHFFFFFCDCEEECCADDEDDDDDDEED	AS:i:-1	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:3	MD:Z:0T99	YT:Z:UU	XS:A:-	NH:i:1
HS21_13853:1:1101:1000:39323#6	147	Schisto_mansoni.Chr_ZW	29878440	50	68M3243N32M	=	29878411	-3372	AATGGATCAAATGAATAACTCAAATCAGCTGTATCCGAATAATAATACTGAATTAGTGATTTTTTATTTTTTCTGCTTGACAAAGAATTTTTGGTATCCG	EDDDDDEEFEEFDDEDDCBAEEEE@EFFFFFFHHJJJJIHIHGJHIHJJJIGHIJJJJJJJJJJJJJJJJJIJJJJIJIIJJJJJJJHHHHHFFFFFCCC	AS:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:100	NM:i:0	XS:A:-	NH:i:1
HS21_13853:1:1101:1000:44900#6	89	Schisto_mansoni.Chr_3	22069811	50	100M	*	0	0	TGTTGGTTTGGCGATGGAGAGAGTGGACTGTGAAGAGGTGTAGTTGTATCGAGAGATGTGAGCGAGCGAGTGTATTTTTGAATTGGAGTGTACTGACTAN	<<B@?<B=B;8?:@A?>>AAB@@;;;;@BBBBBFFECDEEEDCC84=C8F?@BF>IIEEIIGDDDBIIFGEFHFIIIIHIIHEEFG?FC?FFDDBDB:1!	AS:i:-1	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:99G0	YT:Z:UU	NH:i:1	XS:A:+
HS21_13853:1:1101:1000:58791#6	99	Schisto_mansoni.Chr_ZW	51543719	50	19M34N81M	=	51543732	147	NTTCGAGATCCCGTTCCAGGTCTTACTCCAATTCTAGACAAAAGGGGGTCGAACCGTACTATCGTGGTCGCCCTAACGACCGCCGAGATGATTCACGAAG	!1=DDFFFHHHHHIJJIJJJEHIJIJJJJJJJJJJIIDHJJJJIJJJJHGIGHHFFCDDDEFDDCDDDDDDDDDDDDDDDDDDDB@BBDDCEEEDDDDBD	AS:i:-1	XM:i:1	XO:i:0	XG:i:0	MD:Z:0G99	NM:i:1	XS:A:+	NH:i:1
HS21_13853:1:1101:1000:58791#6	147	Schisto_mansoni.Chr_ZW	51543732	50	6M34N94M	=	51543719	-147	TTCCAGGTCTTACTCCAATTCTAGACAAAAGGGGGTCGAACCGTACTATCGTGGTCGCCCTAACGACCGCCGAGATGATTCACGAAGGCGAATAAATAAT	:2BCCAAA>A?:A:C@ACA:>:@:8>DDB>DDB@?<7<?9B?:B@;<9;<?@B>>D?HHEAE;A@@IIIJIJJIJIIGGJIHJJJJJHHGHHFFFFFCCB	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:2	MD:Z:100	YT:Z:UU	XS:A:+	NH:i:1
HS21_13853:1:1101:1001:7799#6	97	Schisto_mansoni.Chr_ZW	17759587	50	99M42N1M	=	17760095	608	NGTCTGACATTTTTTCAGTGCACCTCTCAGAGCATTTTCATTCATATGATAAAATTCTATTTCTTTCTCCAACTCGAAAACCTTTCTGTTTAGGTTTACC	!1=DDFFFHHHHHJJJJJIJJJJJJJJJJJJJJJJJJJIIIJIIIIJJIIJJJJJJIJJIJJJJJJJJJJJJJJJJJHHFFFFFEEEEEEEDDDDDDDDC	AS:i:-1	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:0G99	YT:Z:UU	XS:A:-	NH:i:1
HS21_13853:1:1101:1001:7799#6	145	Schisto_mansoni.Chr_ZW	17760095	50	100M	=	17759587	-608	GAGTATTAAACATCTCTAATAGGTATAAATTCTTGTCTCCAATCTGAGCAACTAAATCTCTGAGTTCATCATTTTCTTTACTGTTATTTTTCAACATTTG	CEEEEEEEFFFFFFFHHHHHHIIJJIJJIJJJJJIIIIIIHGEIIIIJIIJJJJIHFHFHIJJJJJJJJJJJJJIJJJIHGJJJJJJHHHFHFFFFFCCC	AS:i:-6	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:21C78	YT:Z:UU	NH:i:1	XS:A:-
HS21_13853:1:1101:1001:25299#6	99	Schisto_mansoni.Chr_ZW	32388361	50	100M	=	32388529	266	NGTTTTTGTAATTTTGACTGGTGTTACCTCAATTATGGGTAAACGAAGATCTAGATTATTCAATGTTTGTTGTCTTTTTCTTTTGGTAACATAAGGTGAA	!1=DDFFFHHHHHJJJJJJJJIHHIJJJJJJJIJJJJJJFHIIJJIIIHJIIJJIJJJJJHIJJJIJJJHIJHEHHHHFFFFFFEDDEDDDDDDDDCDDC	AS:i:-1	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:0T99	YT:Z:UU	NH:i:1	XS:A:-
HS21_13853:1:1101:1001:25299#6	147	Schisto_mansoni.Chr_ZW	32388529	50	48M2I50M	=	32388361	-266	CATAGGAAAGAATAGAGAAAAATCAAAATATAGATCTCAAAATAAACAATATATATATATATATAGCATAATACTAGAATCTATCAAAAAGTTAAAACTG	FFDECEHHHEHHEIGJHJJIHE@IJJIHGJIHHHFHGIJJIHBJJIGJIGIJJJJJIIJJJJJJJIIIHJIJIIIJJJIHJJJIJJJHHHHHFFFFFB@B	AS:i:-11	XN:i:0	XM:i:0	XO:i:1	XG:i:2	NM:i:2	MD:Z:98	YT:Z:UU	NH:i:1	XS:A:-
HS21_13853:1:1101:1001:48797#6	97	Schisto_mansoni.Chr_1.unplaced.SC_0094	893948	50	33M20753N60M1378N7M	=	916349	22533	NTTTCTTTCGGAATTTTCATATGTATCGATCGCCCTCTGCAACTGTTCCCGTAAGTCTTTTTCCTTCATTTCCAAATCATCTGTAGATAGATACTGTATA	!1=DDFFFHHHHHJJJIJJJJJJJJJJIJJJIIJJIJJJJJIIHJJJFJIIJIIJJIHIIJJIGHHHHHHFDFFFFECEEEEEDEDEEEDDEEDDEDCED	AS:i:-1	XM:i:1	XO:i:0	XG:i:0	MD:Z:0G99	NM:i:1	XS:A:-	NH:i:1

My gtf file is costume made - maybe here's where the problem lies:

Schisto_mansoni.SC_0612	chado	exon	10730	10930	.	+	0	gene_id	"Smp_205760";
Schisto_mansoni.SC_0762	chado	exon	5085	5118	.	-	1	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	5155	5342	.	-	0	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	5632	5868	.	-	0	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	5906	6023	.	-	1	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	6055	6168	.	-	1	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	6208	6285	.	-	1	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	6321	6430	.	-	0	gene_id	"Smp_186560";
Schisto_mansoni.SC_0762	chado	exon	6471	6479	.	-	0	gene_id	"Smp_186560";
Schisto_mansoni.SC_0307	chado	exon	14907	15752	.	+	0	gene_id	"Smp_139870";

Finally, the command I used to run dexseq_counts.py

python ~/bin/dexseq_count.py -p yes -s yes -f bam -r name in_annotation.gtf file.sorted.bam counts.out

I've tried with switching the -s flag but that doesn't seem to make any difference. These same sample I have run before using HTseq_counts and it works a treat.

Any help would be much appreciated.
Thanks,
Anna

[blancha]

Shouldn't the first field of the GTF file correspond to the seqnames?

The third field of the BAM file contains the seqnames Schisto_mansoni.Chr_ZW, Schisto_mansoni.Chr_3, ...
The first field of the GTF file contains the seqnames Schisto_mansoni.SC_0612, Schisto_mansoni.SC_0762, ...

The seqnames do not seem to correspond.
This would seem to be at least one apparent problem.

You can check the proper format for a GTF file here.

http://useast.ensembl.org/info/website/upload/gff.html?redirect=no

The format for a BAM file is defined here.

https://samtools.github.io/hts-specs/SAMv1.pdf

[annaprotasio]

Hi Blancha,

Thanks for your comments.

Agreed, the bam and gtf files do look different. But these are only the ten first lines of the files. Both the bam and the gtf have the information for the mapping and annotation for the whole genome. As far as know, the bam file needs to be sorted by read name and there is no specific order for the gtf.

In addition, I have used the same bam and the same gtf (albeit slightly modified) to run HTseq-counts and it work fine.

Cheers,

Anna

[annaprotasio]

As expected, the problem was in the GFF file.

It was probably me just not reading/interpreting the directions properly, but I suspect this could happen to other people as well.

Be aware that:

1) your GFF (input for dexseq_prepare_annotation.py) MUST contain a transcript_id tag in the last column, as well as a gene_id
2) it is possible that gene_id and transcript_id are expected in that order and not the other way around.

Your input GFF should look like:

Schisto_mansoni.SC_0612	chado	exon	10730	10930	.	+	0	gene_id	"Smp_205760";  transcript_id "Smp_205760.1";
Schisto_mansoni.SC_0762	chado	exon	5085	5118	.	-	1	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	5155	5342	.	-	0	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	5632	5868	.	-	0	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	5906	6023	.	-	1	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	6055	6168	.	-	1	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	6208	6285	.	-	1	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	6321	6430	.	-	0	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0762	chado	exon	6471	6479	.	-	0	gene_id	"Smp_186560"; transcript_id "Smp_186560.1";
Schisto_mansoni.SC_0307	chado	exon	14907	15752	.	+	0	gene_id	"Smp_139870"; transcript_id "Smp_139870.1";