I have few 100 raw fastq files from whole-genome sequencing data and I would like to map these files to a set of genes only (and not whole genome) so as to find SNP's associated with them. Can anyone tell me what snp analysis pipeline is best for targeted variation and not whole genome? I have this in mind: Mapping the reads with bwa and then using varscan Hope to hear for more options
$\begingroup$
$\endgroup$
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against the full reference and then extract the regions you are interested in, e.g. with the region option of samtools view:
samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Check its manual for details.
While there is basically nothign wrong with VarScan it is quite old and not maintained anymore. I suggest you give bcftools a try (check its manual). Other maintained alternatives can be (among others) freebayes, strelka2, VarDict or the GATK from the Broad institute.
$\begingroup$
$\endgroup$
vSNP can be used to target specific regions of a genome. Pipeline is a 2 step process:
- create VCF files (via Freebayes)
- then in your case run all 100 VCF files to create SNP table
A filter file can be supplied to step 2. By filtering all position ranges other than the genes of interest, step 2 will output a SNP table (.xlsx) with only your genes of interest. Also, table will include annotation if GBK is provide.
