Making a Consensus Sequence of a Gene

Question

[User000]

Hello,

I have paired end reads of average 100bp length generated by Illumina. I have a sequence of gene, which is approx. 4000 bp long. I need to produce a consensus sequence of this gene. Can you suggest steps? Is it necessary to do assembly? What I tried to do so far; I mapped my reads directly against the gene using BWA. (Burrows-Wheeler Alignment). Any suggestions are appreciated.

[Ram]

So you have the sequence of the gene, what are you trying to build a consensus of? Also, a consensus needs multiple sequences as input to be called a consensus.

[iraun]

In general I will try to:

1. Mapping the reads against genome
2. Variant calling
3. Variant filtering
4. Extract those variants falling within your gene region (if you know your gene, you can extract the coordinates using bed/gtf/gff3 file)
5. Create consensus (there are different methods and tools, like GATK, vcf-consensus, bcftools consensus... )

[User000]

can I map the reads against the gene of interest, not the whole genome? does it make any sense?

[User000]

My reads come from another organism, and the sequence of gene from another. I followed some steps using BWA, like aligning paired read ends using bwa aln, then combining using bwa sampe, then I used samtools to create a bam file and samtools mpileup to create a consensus.