NGS從比對到群體結(jié)構(gòu)分析軟件及使用方法

醫(yī)學(xué)abeycd 2021-09-24

展開全文

比對分析

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

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

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)計深度信息

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

進行基礎(chǔ)質(zhì)量過濾

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

進行標記群體質(zhì)量過濾

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

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

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

vcftools --vcf vcf --plink --out plink --noweb --ped ped --map map --make-bed --out pro ## 注意plink識別的染色體編號是純數(shù)字，咱們這個是contig的，需要將染色體轉(zhuǎn)換下，后面的數(shù)據(jù)才能正常運行 ## vcftools 在樣品數(shù)量高于 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 # 根據(jù)結(jié)果Ｒ繪制散點圖

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畫圖

軟件文獻

[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