NGS从比对到群体结构分析软件及使用方法 / 开普饭

比对分析

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

排序，因为后续处理都需要按照基因组顺序进行排序

samtools sort bam > out.bam

进行捕获数据的提取和统计

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

去除重复区域

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

统计深度信息

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

建立索引

samtools index bam # 当基因组单条染色体长度大于512M左右的时候，建立索引有问题，后续gatk检测变异会受到影响，建议进行染色体 N 区合适区域进行拆分

变异检测

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按照上面循环下 # 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

NGS从比对到群体结构分析软件及使用方法

相关推荐