==Annotation exercise==
mkdir APOC3 && cd APOC3
vtools init APOC3
vtools import ../data/APOC3.vcf --format ../data/vcf.fmt --build hg19
vtools select 5var --output chr pos ref alt dbNSFP.SIFT_pred dbNSFP.SIFT_score
vtools select 5var --output chr pos ref alt dbNSFP.CADD_raw dbNSFP.CADD_phred
==GATK /IGV exercise==
java -jar GenomeAnalysisTK.jar --help<br data-attributes="%20/" />java -Xmx200m -jar GenomeAnalysisTK.jar -R human_g1k_v37.fa -T UnifiedGenotyper -glm BOTH -I child.bam -I father.bam -I mother.bam -L trio.intervals -o ug.vcf
less -S ug.vcf<br data-attributes="%20/" />java -Xmx200m -jar GenomeAnalysisTK.jar -T VariantEval -R human_g1k_v37.fa -L trio.intervals -D dbsnp.vcf -eval ug.vcf -o ug.vcf.eval<br data-attributes="%20/" />less -S ug.vcf.eval<br data-attributes="%20/" />java -Xmx200m -jar GenomeAnalysisTK.jar -R human_g1k_v37.fa -T SelectVariants -V ug.vcf -mv -mvq 0 -o ug.dnm.vcf -ped trio.ped<br data-attributes="%20/" />java -Xmx500m -jar GenomeAnalysisTK.jar -R human_g1k_v37.fa -T HaplotypeCaller -I child.bam -I father.bam -I mother.bam -minPruning 4 -L 1:199325670-199325672 -mergeVariantsViaLD -o hc.vcf
==PLINKSEQ PSEQ exercise==
Data analysis:
# === Report allele frequencies after 100 generations === #<br data-attributes="%20/" />for (s in s.s) {<br data-attributes="%20/" /> cat("s=") ; cat(s) ; cat(": ")<br data-attributes="%20/" /> freqs = get(paste("freqs.s", s, sep=""))<br data-attributes="%20/" /> cat(freqs[n.gen+1]) ; cat("\n")<br data-attributes="%20/" /> }
# === Graph allele frequency changes === #<br data-attributes="%20/" />pdf("selection_plot.pdf", paper="special", height=4*2, width=4*2, onefile=F)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> for (s in s.s) {<br data-attributes="%20/" /> freqs = get(paste("freqs.s", s, sep=""))<br data-attributes="%20/" /> lines(0:n.gen,freqs, col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" />dev.off()
==<br data-attributes="%20class=%22Apple-interchange-newline%22">Regression exercise==
In R:
<br data-attributes="%20/" />spower show exercise.SEQPowerDB LOGIT method power title --condition "where power between 0.25 and 0.95"<br data-attributes="%20/" />
<br data-attributes="%20/" />for i in 1 1.5 2 2.5 3 3.5 4; do<br data-attributes="%20/" />spower LOGIT Kryukov2009European1800.sfs --sample_size 1000 --OR_rare_detrimental $i --method "CFisher --name CMC$i" --title FixedOR$i -r 100 -j 4 -l 1 -o exercise2<br data-attributes="%20/" />done
==VAT exercise==
vtools -h
vtools remove variants to_remove -v0
vtools show tables
vtools remove genotypes "DP_geno<10" -v0 <br /> vtools select variant "mut_type like 'non%' or mut_type like 'stop%' or region_type='splicing'" -t v_funct <br /> vtools show tables <br /> vtools show samples --limit 5 <br /> vtools select variant --samples "RACE=1" -t CEU <br /> mkdir -p ceu
cd ceu
vtools init ceu --parent ../ --variants CEU --samples "RACE=1" --build hg19
vtools show project
vtools select variant "CEU_mafGD10>=0.05" -t common_ceu
vtools select v_funct "CEU_mafGD10<0.01" -t rare_ceu <br /> vtools use refGene <br /> vtools show annotation refGene <br /> vtools associate -h <br /> vtools show tests <br /> vtools show test LinRegBurden
vtools associate common_ceu BMI --covariate SEX -m "LinRegBurden --alternative 2" -j1 --to_db EA_CV > EA_CV.asso.res
grep -i error *.log
vtools show fields
vtools associate rare_ceu BMI --covariate SEX -m "LinRegBurden --alternative 2" -g refGene.name2 -j1 --to_db EA_RV > EA_RV.asso.res
grep -i error *.log | tail -2210
less EA_RV.asso.res
sort -g -k6 EA_RV.asso.res | head
vtools associate rare_ceu BMI --covariate SEX -m "VariableThresholdsQt --alternative 2 -p 100000 --adaptive 0.0005" -g refGene.name2 -j1 --to_db EA_RV > EA_RV_VT.asso.res
grep -i error *.log | tail -2210
less EA_RV_VT.asso.res
sort -g -k6 EA_RV_VT.asso.res | head
vtools select rare_ceu "refGene.name2='ABCC1'" -o chr pos ref alt CEU_mafGD10 numGD10 mut_type --header
vtools_report plot_association qq -o QQRV -b --label_top 2 -f 6 < EA_RV.asso.res
vtools_report plot_association manhattan -o MHRV -b --label_top 5 --color Dark2 --chrom_prefix None -f 6 < EA_RV.asso.res <br /> vtools associate rare_ceu BMI --covariate SEX KING_MDS1 KING_MDS2 -m "LinRegBurden --name RVMDS2 --alternative 2" -g refGene.name2 -j1 --to_db EA_RV > EA_RV_MDS2.asso.res vtools_report plot_association qq -o QQRV_MDS2 -b --label_top 2 -f 6 < EA_RV_MDS2.asso.res <br /> cd .. <br /> vtools select variant --samples "RACE=0" -t YRI<br /> mkdir -p yri<br />; cd yri
vtools init yri --parent ../ --variants YRI --samples "RACE=0" --build hg19
vtools select variant "YRI_mafGD10>=0.05" -t common_yri
vtools select v_funct "YRI_mafGD10<0.01" -t rare_yri <br /> vtools use refGene vtools associate common_yri BMI --covariate SEX -m "LinRegBurden --alternative 2" -j1 --to_db YA_CV > YA_CV.asso.res
vtools associate rare_yri BMI --covariate SEX -m "LinRegBurden --alternative 2" -g refGene.name2 -j1 --to_db YA_RV > YA_RV.asso.res
vtools associate rare_yri BMI --covariate SEX -m "VariableThresholdsQt --alternative 2 -p 100000 --adaptive 0.0005" -g refGene.name2 -j1 --to_db YA_RV > YA_RV_VT.asso.res
cd ..
vtools_report meta_analysis ceu/EA_RV_VT.asso.res yri/YA_RV_VT.asso.res --beta 5 --pval 6 --se 7 -n 2 --link 1 > META_RV_VT.asso.res
cut -f1,3 META_RV_VT.asso.res | head==VMT==