比對(duì)分析

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進(jìn)行基因組建index

排序，因?yàn)楹罄m(xù)處理都需要按照基因組順序進(jìn)行排序

samtools sort bam > out.bam

進(jìn)行捕獲數(shù)據(jù)的提取和統(tǒng)計(jì)

samtools view -@ 6 -L  bed -b -o bam input.bam samtools flagstat bam > flagstat

去除重復(fù)區(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)計(jì)深度信息

samtools depth -q 20 -Q 20 -l 60 -d 5000 bam > depth

建立索引

samtools index bam # 當(dāng)基因組單條染色體長(zhǎng)度大于512M左右的時(shí)候，建立索引有問題，后續(xù)gatk檢測(cè)變異會(huì)受到影響，建議進(jìn)行染色體 N 區(qū)合適區(qū)域進(jìn)行拆分

變異檢測(cè)

gatk進(jìn)行單個(gè)樣品gvcf檢測(cè)

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

進(jìn)行基礎(chǔ)質(zhì)量過(guò)濾

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

進(jìn)行標(biāo)記群體質(zhì)量過(guò)濾

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   ## 建議有條件的測(cè)序 10X，個(gè)體深度>=7,miss 0.1；maf 0.05，不足的話，4X  

　群體結(jié)構(gòu)分析

數(shù)據(jù)格式轉(zhuǎn)化

vcftools --vcf vcf --plink --out plink  --noweb --ped ped --map map --make-bed --out pro ## 注意plink識(shí)別的染色體編號(hào)是純數(shù)字，咱們這個(gè)是contig的，需要將染色體轉(zhuǎn)換下，后面的數(shù)據(jù)才能正常運(yùn)行 ## vcftools 在樣品數(shù)量高于 1008 左右的時(shí)候會(huì)有問題，可以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 # 根據(jù)結(jié)果Ｒ繪制散點(diǎn)圖

TREE 分析

# 將VCF|PED文件轉(zhuǎn)換為fa格式 treebest nj -b 1000  fa > nj.tree ggtree 或者其他在線軟件，離線軟件畫圖

STRUCTURE分析

admixture  --cv 1 > 1.log ## 1-10按照上面循環(huán)下 # R畫圖

軟件文獻(xiàn)

[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.

 作者：點(diǎn)滴生信 https://www.bilibili.com/read/cv6573709/ 出處：bilibili