Post date: May 13, 2014 9:54:40 PM
#Need to index sam version of genome:
samtools faidx Esosorum_40mil_qc.mmp84.fasta
#See jobsubVariantCall.sh for variant calling script. Changed path and names. Outfiles will be varEsosorum.bcf, varEsosorum.vcf (this one will be readable). To run:
qsub jobsubVariantCall.sh
#254919.torque
#See wrap_qsub_rc_theta.pl for getting theta and h for each of 9 sosorum subpopulations: C-EL, C-F1, C-F2, C-F3, W-EL, W-CS, W-PO, W-SG, W-UP. We may combine wild pops into one in the future. To run:
perl wrap_qsub_rc_theta.pl
cd /labs/evolution/projects/sosorum/
samtools mpileup -C 50 -d 50 -I -g -q 10 -Q 15 -f Esosorum_40mil_qc_mmp84.fasta alnE-BS-W-UP*sorted.bam > W-UPsites.bcf
bcftools view -cGP cond2 W-UPsites.bcf > /dev/null 2> W-UPsites.1.afs
bcftools view -cGP W-UPsites.1.afs W-UPsites.bcf > /dev/null 2> W-UPsites.2.afs
bcftools view -cGP W-UPsites.2.afs W-UPsites.bcf > /dev/null 2> W-UPsites.3.afs
bcftools view -cGP W-UPsites.3.afs W-UPsites.bcf > /dev/null 2> W-UPsites.4.afs
bcftools view -cGP W-UPsites.4.afs W-UPsites.bcf > /dev/null 2> W-UPsites.5.afs
bcftools view -cGP W-UPsites.5.afs W-UPsites.bcf > /dev/null 2> W-UPsites.6.afs
bcftools view -cGP W-UPsites.6.afs W-UPsites.bcf > /dev/null 2> W-UPsites.7.afs
bcftools view -cGP W-UPsites.7.afs W-UPsites.bcf > /dev/null 2> W-UPsites.8.afs
bcftools view -cGP W-UPsites.8.afs W-UPsites.bcf > /dev/null 2> W-UPsites.9.afs
bcftools view -cGP W-UPsites.9.afs W-UPsites.bcf > /dev/null 2> W-UPsites.10.afs
bcftools view -cGP W-UPsites.10.afs W-UPsites.bcf > /dev/null 2> W-UPsites.11.afs
bcftools view -cGP W-UPsites.11.afs W-UPsites.bcf > /dev/null 2> W-UPsites.12.afs
bcftools view -cGP W-UPsites.12.afs W-UPsites.bcf > /dev/null 2> W-UPsites.13.afs
bcftools view -cGP W-UPsites.13.afs W-UPsites.bcf > /dev/null 2> W-UPsites.14.afs
bcftools view -cGP W-UPsites.14.afs W-UPsites.bcf > /dev/null 2> W-UPsites.15.afs
bcftools view -cGP W-UPsites.15.afs W-UPsites.bcf > /dev/null 2> W-UPsites.16.afs
bcftools view -cGP W-UPsites.16.afs W-UPsites.bcf > /dev/null 2> W-UPsites.17.afs
bcftools view -cGP W-UPsites.17.afs W-UPsites.bcf > /dev/null 2> W-UPsites.18.afs
bcftools view -cGP W-UPsites.18.afs W-UPsites.bcf > /dev/null 2> W-UPsites.19.afs
bcftools view -cGP W-UPsites.19.afs W-UPsites.bcf > /dev/null 2> W-UPsites.20.afs
#Will use a file that looks like W-UPsite20.vcf for each subpop.
#255163-255171.torque
#May 15: The 20 runs of h/ theta did not look like they converged by the end, so we decided to run 10 more. We modified the wrap_qsub_rc_theta.pl script to do this.
#256701-256709.torque
#Filtered the .vcf file for quality control, see: vcfFilter.pl
#mincoverage = 774 (this is 3 x 258 to get 3x coverage per individual per locus)
#minAltRds = 10
#maxRevOrient = 0
#changed the word scaffold to Contig
perl vcfFilter.pl varEsososrum.vcf
#retained 24507 variable loci
perl getCoverage.pl filtered_varEsosorum.vcf
less coverage.txt #this is the coverage at each locus
#read into R to get average coverage
>mean (co[,1])/258
#average coverage is 5.99 x