Creating a mRNA GTF file from fasta for HTSeq ?

[CEPHSeq]

Hi,

I have a sam file (bwa) from paired end RNASeq short reads that was aligned to a CDS fasta file. I need to use HTSeq.counts on this sam file, so I need the corresponding GTF file.

I thought I would not be too had to generate a basic GTF file from the original RNA fasta file, but HTSeq does not recognize the ID of any RNA sequence in the SAM file :

Warning: Skipping read 'XXX0654:58235#TGACCA', because chromosome 'gi|155030243|ref|NM_017599.3|', to which it has been aligned, did not appear in the GFF file.

However, this sequence id is present in the GTF file :
gi|155030243|ref|NM_017599.3| ref CDS 1 4580 . + . gene_id "VEZT"; transcript_id "NM_017599.3";
and in the sam header too :
@SQ SN:gi|155030243|ref|NM_017599.3| LN:4580

Where is the mismatch ??

Many thanks for your help.

Emmanuel.

[CEPHSeq]

Sorry, I answer myself : bad reference version. My gtf is OK but the rna.fasta is from another release, so the NM_XXX accession version was different by the number
after the dot, hence the mismatch.

Emmanuel.

[sgcsd]

Hi,
Can you provide the info on how you created GTF file from fasta sequence file..

Thank you

[kurban910]

i wanna repeat the same question, because i am working with the transcripome data of non model organism, but during downstream analysis i am often required GFT files of the fasta file i work with, since this insect does not have a genomic information available, so i could not get the transcripts coordinate data. am thinking whether or not i can creat GFT file, that could enable me to count reads with HTseq-count and edgeR DEG analysis.
if it could be done, how?

i do not know my question is ridiculous or not, but i want you guys suggestion.

[dpryan]

When you align to the transcriptome like this you don't use htseq-count. Instead, filter out secondary alignments and anything else you want (e.g., remove alignments with very low mapping quality), index the results, and use "samtools idxstats" to get the counts. Alternatively, you can use RSEM or eXpress (or one of the many equivalents) to get estimated counts.

[kurban910]

hi @dpryan,
thanks for the reply. but what do u mean by "index the results"? do u mean index the sam file i got after aligning the reads to trinity.fasta by using bowtie2?

does that mean i could not analyze my data by using edgeR, DESeq, and baySeq ?

[dpryan]

You'll have to convert it to BAM and sort it, but yes.

You'll get unique integer counts with the idxstats method, so edgeR/DESeq2/etc. will work fine. For the RSEM/eXpress/etc. route, edgeR will work but DESeq2 will not (no clue about baySeq).

[kurban910]

thank you @Devon Ryan, really

[kurban910]

hi @Devon,
i have finished differential the expression analysis with edgeR with the help of its user guide, then searched some stuff about how to export data from R. but i still do not know how to export the edgeR analysis results, am using Ubuntu . could you give me some basic tips there?

[dpryan]

The "write.table()" function is probably the most convenient method. I assume you have some sort of data frame that you'd like written to a text file so you/others can easily use it (e.g., in Excel, as supplemental data, or again in R).

[dpryan]

Either read it into R yourself and then deal with the columns however you'd like or use awk to create a new text file of each. You probably want to either just count read #1 or filter out singleton alignments, which can be done with samtools.

[kurban910]

hey @dpryan,
it seems like i did not make myself clear, sorry about it.
i used edgeR got this:

et <- exactTest(y, dispersion=bcv^2)
> summary(de <- decideTestsDGE(et))
   [,1] 
-1   273
0  27700
1    877
>

is there any function that i could use to get the 273 down-regulated and 877 up-regulated transcripts?