count file having zero

Question

[BioBaby]

I have some non-coding region sequence around 5000 and some tissue specific reads data. I want to do differential gene expression analysis. So I want to check which tissue specific reads having these non coding sequences for that I aligned them (where I created database for those non coding sequences) by tophat2. Then I need to have count file from these bam file, used featureCount but it gave zero only.I also used bedtool also it also not worked. I suspect because of those non-coding seq it can't generate count file...may be I don't know. I convert bam to sam file then it looks like this:

> @HD     VN:1.0  SO:coordinate
@SQ     SN:0    LN:309
@SQ     SN:0    LN:381
@SQ     SN:1    LN:1085
@SQ     SN:1    LN:353
@SQ     SN:10   LN:239
@SQ     SN:10   LN:205
@SQ     SN:10008        LN:659
@SQ     SN:1001 LN:291

Can you people tell me what should I do?

[Benn]

I changed "tool" to "question", since tool is meant for promotion of new tools. Also use the 101 button for code, it reads easier for your readers.

[ATpoint]

Your post lacks any details (code) on what you did and which files you used). Please add this. Brief Reminder On How To Ask A Good Question

[Devon Ryan]

You first need to figure out why you have duplicate entries with different lengths in your BAM/SAM files. For example, you an a contig 0 and length 309 and another contig 0 with length 381. That's going to break most tools. Further, make sure that your annotations has contig names like 0, 1, 10008 and so on, otherwise featureCounts will always return 0.

[BioBaby]

Thanks Ryan...I removed duplicate. Here I should mention that I am working novel lincRNA so my annotation file does not have that information.Then how to calculate count? any other way to do it. It may look silly but do not have idea.

[Devon Ryan]

Add it to your annotation then.

[BioBaby]

here I am giving my gtf file which looks like:

SL3.0ch00       pfurio  peak    0       16479   .       +       .       peak_id "SL3.0ch00_0_16479";
SL3.0ch00       pfurio  peak    17940   328351  .       +       .       peak_id "SL3.0ch00_17940_328351";
SL3.0ch00       pfurio  peak    334459  548343  .       +       .       peak_id "SL3.0ch00_334459_548343";

from bam file samid generated by following command

samtools idxstats a_sorted.bam > samid.txt

samid.txt looks like this:

SL3.0ch00:1230877-1231088       211     0       0
SL3.0ch00:1147815-1150318       2503    0       0
SL3.0ch00:1147485-1147759       274     0       0
SL3.0ch00:1150608-1151254       646     4       0

then give this command

featureCounts -a x.gtf -F GTF -o y.bed a_sorted.bam

ERROR: no features were loaded in format GTF

and I used bedtool also still it gave zero count.

[Devon Ryan]

That "GTF" file isn't in GTF format, you need to fix that. Also, the chromosome names don't match between the two files.

[BioBaby]

I created GTF file from this command:

awk -v OFS="\t" '{print $1, "myIntergenic","intergenic", $2, $3, ".", "+", ".", sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"" ";", $4, $1, $2, $3 )}' l.bed > o.gtf

so it looks like this:

SL3.0ch00   myIntergenic    intergenic  1230877 1231088 .   +   .   gene_id "xxx"; transcript_id "SL3.0ch00:1230877-1231088";
SL3.0ch00   myIntergenic    intergenic  1147815 1150318 .   +   .   gene_id "xxx"; transcript_id "SL3.0ch00:1147815-1150318";
SL3.0ch00   myIntergenic    intergenic  1147485 1147759 .   +   .   gene_id "xxx"; transcript_id "SL3.0ch00:1147485-1147759";

but how to change bam file ids? in gtf file we can't give "SL3.0ch00:1230877-1231088"...so what should I do?

[Devon Ryan]

It looks like your reference fasta file has the contigs split apart into chunks. Honestly I'd try to fix that and realign the data. Then you'd have SL3.0ch00 and such and the alignments would probably be more accurate.

[BioBaby]

thanks for giving your time Ryan...one question I have..can I create count file from idxstats output? it will be OK for further differential gene expression analysis?

[Devon Ryan]

In general that won't produce correct results.