RNA read alingment to draft genome

[vinaydu]

Dear all,

I am working on a plant system. I have sequenced the transcriptome of the two species of the same genus. Now a draft genome has been released recently for a third species. I want to find out SNPs in case of these three species. Here is my approach that I have taken.

1. Align rna-seq paired reads using STAR aligner to the draft genome. (>85% alignment I am getting after playing with various options).

2. Convert the alignment into BAM.

3. Sort the alignment.

4. SNP call using smatools mpileup and bcftools.

Here I am getting an error Incorrect number of fields!!!!!

Could it has anything to do with the extra information in the header of the reference draft genome???


Please help.....

[westerman]

Step 2 is probably your problem. How are you converting to BAM? Not using STAR I am unsure if it can output BAM natively. If so then that is your best bet.

[vinaydu]

@westerman

Thanks for replying.

Yes....I am cnverting sam to bam by using samtools.

Tough in the current version of STAR direct output in BAM format is also possible.

So what you suggest???

Shall I not be using smatools to convert to BAM ?
Should I get the input directly in BAM format?

Though I have again started this conversion using samtools. I shall be posting the results once it finished.

[alexdobin]

You can now do both conversion to BAM and BAM sorting directly from STAR using:
--outSAMtype BAM SortedByCoordinate