比對分析
bwa
bwa mem -t 6 -k 32 -M -R "@RG\tID:saample\tLB:sample\tSM:sample" fa fq_R1.fq.gz fq_R2.fq.gz |samtools view -b -S ->sample.bam # bwa index進行基因組建index
排序，因為后續(xù)處理都需要按照基因組順序進行排序
samtools sort bam > out.bam
進行捕獲數據的提取和統(tǒng)計
samtools view -@ 6 -L bed -b -o bam input.bam samtools flagstat bam > flagstat
去除重復區(qū)域
java -jar picard-1.119/MarkDuplicates.jar I=bam O=rmdup.bam M=.metrics \ VALIDATION_STRINGENCY=SILENT CREATE_INDEX=true REMOVE_DUPLICATES=true TMP_DIR=tmp MAX_FILE_HANDLES=1000
統(tǒng)計深度信息
samtools depth -q 20 -Q 20 -l 60 -d 5000 bam > depth
建立索引
samtools index bam # 當基因組單條染色體長度大于512M左右的時候，建立索引有問題，后續(xù)gatk檢測變異會受到影響，建議進行染色體 N 區(qū)合適區(qū)域進行拆分
變異檢測
gatk進行單個樣品gvcf檢測
gatk-4.1.2.0/gatk HaplotypeCaller -R ref --emit-ref-confidence GVCF -I rmdup.bam -O g.vcf # 使用前將基因組建好索引 gatk CreateSequenceDictionary -R fa
合并gvcf
gatk-4.1.2.0/gatk CombineGVCFs -R ref -O merge.g.vcf -V gvcf ...
call 基因型
gatk GenotypeGVCFs -R ref -V merge.g.vcf -O raw.vcf
拆分snp，indel
gatk SplitVcfs --INPUT=raw.vcf --INDEL_OUTPUT=raw.indel.vcf --SNP_OUTPUT=raw.SNP.vcf --STRICT=false
進行基礎質量過濾
gatk-4.1.2.0/gatk VariantFiltration -V raw.SNP.vcf \ -filter "QD < 2.0" --filter-name "QD2" \ -filter "QUAL < 30.0" \ --filter-name "QUAL30" \ -filter "SOR > 3.0" \ --filter-name "SOR3" \ -filter "FS > 60.0" \ --filter-name "FS60" \ -filter "MQ < 40.0" \ --filter-name "MQ40" \ -filter "MQRankSum < -12.5" \ --filter-name "MQRankSum-12.5" \ -filter "ReadPosRankSum < -8.0" \ --filter-name "ReadPosRankSum-8" -O raw.SNP.filter.QD2.QUAL30.SOR3.FS60.MQ40.MQRankSum-12.5.ReadPosRankSum-8.vcf gatk-4.1.2.0/gatk VariantFiltration \ -V raw.indel.vcf \ -filter "QD < 2.0" \ --filter-name "QD2" \ -filter "QUAL < 30.0" \ --filter-name "QUAL30" \ -filter "FS > 200.0" \ --filter-name "FS200" \ -filter "ReadPosRankSum < -20.0" \ --filter-name "ReadPosRankSum-20" \ -O raw.indel.QD2.QUAL30.FS200.ReadPosRankSum-20.vcf
進行標記群體質量過濾
vcftools --minDP 4 --maxDP 100 --minGQ 10 --minQ 30 --min-meanDP 3 \ --out aw.SNP.filter.QD2.QUAL30.SOR3.FS60.MQ40.MQRankSum-12.5.ReadPosRankSum-8.minDP4.maxDP100.minGQ10.minQ30.min-meanDP3.miss0.2.maf0.01.vcf \ --vcf raw.SNP.filter.QD2.QUAL30.SOR3.FS60.MQ40.MQRankSum-12.5.ReadPosRankSum-8.vcf \ --recode --recode-INFO-all --max-missing 0.8 --maf 0.01 ## 建議有條件的測序 10X，個體深度>=7,miss 0.1；maf 0.05，不足的話，4X
群體結構分析
數據格式轉化
vcftools --vcf vcf --plink --out plink --noweb --ped ped --map map --make-bed --out pro ## 注意plink識別的染色體編號是純數字，咱們這個是contig的，需要將染色體轉換下，后面的數據才能正常運行 ## vcftools 在樣品數量高于 1008 左右的時候會有問題，可以1000 樣品一組生成 ped 和 map，然后直接 將 ped cat合并在一起即可，理論上，map 是一致的
PCA分析
gcta64 --bfile pro --make-grm --thread-num 4 --out pro gcta64 --grm pro --pca 10 --thread-num 4 --out pro # 根據結果Ｒ繪制散點圖
TREE 分析
# 將VCF|PED文件轉換為fa格式 treebest nj -b 1000 fa > nj.tree ggtree 或者其他在線軟件，離線軟件畫圖
STRUCTURE分析
admixture --cv 1 > 1.log ## 1-10按照上面循環(huán)下 # R畫圖
軟件文獻
[1] Li, H. and R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009. 25(14): p. 1754-60.
[2] Li, H., et al., The Sequence Alignment/Map format and SAMtools.Bioinformatics, 2009. 25(16): p. 2078-9.
[3] Danecek, P., et al., The variant call format and VCFtools. Bioinformatics, 2011. 27(15): p. 2156-8.
[4] Purcell, S., et al., PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet, 2007. 81(3): p. 559-75.
[5] Yang, J., et al., GCTA: A Tool for Genome-wide Complex Trait Analysis. The American Journal of Human Genetics, 2011. 88(1): p. 76-82.
[6] Alexander, D.H., J. Novembre, and K. Lange, Fast model-based estimation of ancestry in unrelated individuals. Genome Res, 2009. 19(9): p. 1655-64.
[7] Vilella, A.J., et al., EnsemblCompara GeneTrees: Complete, duplication-aware phylogenetic trees in vertebrates. Genome Res, 2009. 19(2): p. 327-35.
作者：點滴生信 https://www.bilibili.com/read/cv6573709/ 出處：bilibili
