Post date: May 15, 2014 4:49:28 PM
reuse SAMtools
samtools faidx Enana_40mil_qc.mmp84.fasta
#See jobsubVariantCall.sh for variant calling script. Changed path and names. Outfiles will be varEnana.bcf, varEnana.vcf (this one will be readable). To run:
qsub jobsubVariantCall.sh
#256717.torque
#See wrap_qsub_rc_theta.pl for getting theta and h for each of 7 nana subpopulations: C-1, C-2, C-3, C-4, W-BD, W-DS, W-HO. We may combine wild pops into one in the future. To run:
perl wrap_qsub_rc_theta.pl
#256734-256740.torque
#Filtered the .vcf file for quality control, see: vcfFilter.pl
#mincoverage = 585 (this is 3 x 195 to get 3x coverage per individual per locus)
#minAltRds = 10
#maxRevOrient = 0
perl vcfFilter.pl varEnana.vcf
#retained 456 variable loci
perl getCoverage.pl filtered_varnana.vcf
less coverage.txt #this is the coverage at each locus
#read into R to get average coverage
>mean (co[,1])/195
#average coverage is 10.22 x