一、数据准备
1、基因组：
dh_hic.fasta
2、重测序数据：
illumina双端测序，经过fastp质控和过滤之后：
BBL1_clean_1.fq.gz、BBL1_clean_2.fq.gz
BBL2_clean_1.fq.gz、BBL2_clean_2.fq.gz
... ...
#共计129个个体，258个文件
3、基因组注释文件：
EVM.dh.gff
4、目录结构：
/re-sequence/
├── /1.variant_calling #鉴定SNP
├── /2.vairant_annotation #变异注释
├── /3.polulation_genetics #群体遗传学分析
└── /software
二、SNP calling
/re-sequence/
├── /1.variant_calling
├── /1.ref
├── /2.map_resequenceDATA_to_ref
├── /3.SNP_and_InDel
1、在 /1.ref/ 文件夹中，对参考基因组构建索引
samtools faidx dh_hic.fasta
bwa index dh_hic.fasta
java -jar /home/lx_sky6/yt/soft/picard.jar CreateSequenceDictionary R=dh_hic.fasta
2、在 /2.map_resequenceDATA_to_ref/ 文件夹中，使用bwa软件将重测序数据分别比对到参考基因组
(1) 命令如下，举例其中一条：
bwa mem -t 30 -R '@RG\tID:BBL1\tSM:BBL1\tPL:illumina' \ #read group header line
/home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \ #参考基因组
/public/lx_sky/ayu/resequencing/data/clean_data/BBL1_clean_1.fq.gz \ #双端测序的文件1
/public/lx_sky/ayu/resequencing/data/clean_data/BBL1_clean_2.fq.gz \ #双端测序的文件2
| /home/lx_sky6/software/miniconda3/envs/rna-seq/bin/samtools sort -@ 30 -m 100G -o BBL1.sort.bam - #使用samtools对比对生成的.sam文件进行排序，并转换为.bam格式
运行此步骤后，129个个体生成了129个 XXX.sort.bam 文件
（2） 使用 picard.jar 软件包中的 MarkDuplicates标记重复序列
命令如下，举例其中一条：
java -Xmx250g -XX:ParallelGCThreads=63 -jar /home/lx_sky6/yt/soft/picard.jar MarkDuplicates \
I=BBL1.sort.bam \ #输入上一步生成的XXX.sort.bam
O=BBL1.sort.rmdup.bam \ #输出文件名
CREATE_INDEX=true REMOVE_DUPLICATES=true \ #创建索引并删除重复序列信息
M=BBL1.marked_dup_metrics.txt
运行此步骤后，129个个体生成了129个 XXX.sort.rmdup.bam 文件
（3）统计比对率，命令举例：
samtools flagstat BBL1.sort.bam > BBL1.sort.bam.flagstat
3、在 /3.SNP_and_InDel/ 文件夹中，进行SNP和InDel的鉴定、筛选
（1）使用 gatk 软件包的 HaplotypeCaller 模块鉴定SNP、InDel
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" HaplotypeCaller \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \ #参考基因组
-I /home/lx_sky6/liuying/re-sequence/1.variant_calling/2.map_resequenceDATA_to_ref/BBL1.sort.rmdup.bam \ #上一步生成的XXX.sort.rmdup.bam文件
-ERC GVCF \ #生成gvcf格式文件
-O BBL1.g.vcf \ #每个个体的gvcf文件名
1>BBL1.HC.log 2>&1 #日志文件
运行此步骤后，129个个体生成了129个 XXX.g.vcf 文件
(2) 使用 gatk 软件包的 CombineGVCFs 模块合并gvcf文件
ls ./*.g.vcf > gvcf.list #获取129个 .g.vcf 文件列表
#运行以下命令
gatk --java-options "-Xmx150g -Djava.io.tmpdir=./tmp" CombineGVCFs \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V gvcf.list \
-O dh.merge.g.vcf
运行此步骤后，生成了1个 dh.merge.g.vcf文件，存放着129个个体的SNP、InDel信息
（3） 使用 gatk 软件包的 GenotypeGVCFs 模块，对dh.merge.g.vcf文件进行基因分型
gatk --java-options "-Xmx150g -Djava.io.tmpdir=./tmp" GenotypeGVCFs \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
--variant dh.merge.g.vcf \
-O dh.merge_raw.vcf
运行此步骤后，生成了1个 dh.merge_raw.vcf文件
（4） 挑选出SNP、并进行初步硬过滤
#step1：从dh.merge_raw.vcf文件挑选出SNPs
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" SelectVariants \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V dh.merge_raw.vcf \
--select-type SNP \
-O dh.raw.snp.vcf
#step2：设置过滤条件，给 满足/不满足 条件的SNPs添加标签
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" VariantFiltration \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta -V dh.raw.snp.vcf \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name 'SNP_filter' \
-O dh.filter.snp.vcf
#step3：只保留满足过滤条件的SNPs位点信息
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" SelectVariants \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V dh.filter.snp.vcf \
--exclude-filtered \
-O dh.filtered.snp.vcf
运行此步骤后，满足过滤条件的SNPs被储存到 dh.filtered.snp.vcf 文件中
（5） 挑选出InDel、并进行初步筛选、过滤
#step1：从dh.merge_raw.vcf文件挑选出InDels
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" SelectVariants \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V dh.merge_raw.vcf \
--select-type INDEL \
-O dh.raw.indel.vcf
#step2：设置过滤条件，给 满足/不满足 条件的InDels添加标签
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" VariantFiltration \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V dh.raw.indel.vcf \
--filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name 'INDEL_filter' \
-O dh.filter.indel.vcf
#step3：只保留满足过滤条件的Indels位点信息
gatk --java-options "-Xmx100g -Djava.io.tmpdir=./tmp" SelectVariants \
-R /home/lx_sky6/liuying/re-sequence/1.variant_calling/1.ref/dh_hic.fasta \
-V dh.filter.indel.vcf \
--exclude-filtered \
-O dh.filtered.indel.vcf
运行此步骤后，满足过滤条件的InDels被储存到 dh.filtered.snp.vcf 文件中
三、变异检测结果注释
/re-sequence/
├── /1.variant_calling
├── /2.vairant_annotation #变异注释
├── /3.polulation_genetics
1、待添加
四、群体遗传分析
/re-sequence/
├── /1.variant_calling
├── /2.vairant_annotation
├── /3.polulation_genetics #群体遗传学分析
├── /01.filter
├── /02.vcf_to_phylip
├── /1.Tree
├── /2.PCA
├── /3.Structure
├── /4.LDdecay
└── /5.positive_selection
1、对SNP的vcf文件进行进一步过滤
（1）在 /01.filter/ 文件中：
plink --vcf dh.filtered.snp.vcf \ #
--geno 0.1 \ #允许某个SNP在所有样本中的缺失率最大为10%，如：100个体中，有11个个体在该位点信息缺失，则舍去此SNP位点
--maf 0.05 \ #该位点的等位基因频率，最小值为5%。
#如一个SNP位点存在A和C两种，那么A和C的最小基因频率都应该大于5%。
#之所以用这个过滤标准，是因为MAF如果非常小，比如低于0.05，
#那么意味着大部分位点都是相同的基因型，这些位点贡献的信息非常少，增加假阳性
--out dh.SNP.missing_maf \
--recode vcf-iid \
--allow-extra-chr \
--set-missing-var-ids @:# \
--keep-allele-order
#另外：还可用vcftools过滤掉非【双等位基因SNP】，如下：
vcftools --gzvcf combined200.vcf.gz --min-alleles 2 --max-alleles 2
#适用于二倍体物种？
得到的 dh.SNP.missing_maf 可用于：
———正选择分析（Fst、Pi、tajimaD、ROD-Fst、xpclr）
———LD衰减分析
（注意 dh.SNP.missing_maf 文件的适用范围，区分于 dh.LDfilter.vcf ）
（2） 在 dh.SNP.missing_maf 的基础上，进行LD位点过滤
plink --vcf dh.SNP.missing_maf \
--indep-pairwise 50 10 0.2 \
--out tmp.ld \
--allow-extra-chr \
--set-missing-var-ids @:# &
#然后接着运行：
plink --vcf dh.SNP.missing_maf \
--make-bed \
--extract tmp.ld.prune.in \
--out dh.LDfilter \
--recode vcf-iid \
--keep-allele-order \
--allow-extra-chr \
--set-missing-var-ids @:# &
得到过滤掉LD位点后的 dh.LDfilter.vcf 文件，可用于：
———进化树分析
———PCA
———Structure分析
（注意 dh.LDfilter.vcf 文件的适用范围，区分于 dh.SNP.missing_maf ）
