###############################################################
####### #######
####### 1 reads比對 #######
####### 20220108 #######
###############################################################

基因組文件：reference.fasta
測序數(shù)據(jù)文件：example_R1.fq.gz example_R2.fq.gz

1 為參考基因組構(gòu)建索引

samtools faidx reference.fasta
bwa index reference.fasta

2 將測序數(shù)據(jù)文件比對到參考基因組

-t 設(shè)置線程數(shù),-R 添加reads的標(biāo)簽信息,view -Sb 將結(jié)果文件由sam轉(zhuǎn)為bam，減少存儲。

bwa mem -t 2 \ 
-R '@RG\tID:example\tPL:illumina\tSM:example' \ 
reference.fasta \ 
./example_R1.fq.gz ./example_R2.fq.gz \ 
| samtools view -Sb - > example.bam  

3 比對結(jié)果文件進(jìn)行排序

-@ 設(shè)置線程數(shù)，-m 每個線程可使用的最大內(nèi)存，-o 指定輸出文件名稱

samtools sort -@ 2 -m 10G -o example.sort.bam example.bam

4 去除比對結(jié)果中的PCR重復(fù)

samtools rmdup -S example.sort.bam example.sort.rmdup.bam

5 過濾比對質(zhì)量少于20的reads

samtools view example.sort.rmdup.bam -q 20 -O BAM -o example.sort.rmdup.filter.bam

###############################################################
####### #######
####### 2 變異檢測 #######
####### 20220108 #######
###############################################################

基因組文件：reference.fasta
比對結(jié)果文件：example.bam

注：GATK3中需要單獨(dú)對bam文件進(jìn)行局部重比對，以提高Indel附近區(qū)域比對結(jié)果的準(zhǔn)確性，但在GATK4中，

局部重比對的步驟已經(jīng)自動包含在HaplotypeCaller模塊中，無需單獨(dú)進(jìn)行。

1 分別對每個樣品進(jìn)行變異檢測

用到的GATK模塊為HaplotypeCaller，--java-options 設(shè)置JAVA的參數(shù) -XX:+UseSerialGC 單線程運(yùn)行，-ERC 輸出GVCF

gatk --java-options "-XX:+UseSerialGC"
HaplotypeCaller -R reference.fasta
-ERC GVCF
-I example.bam
-O example.g.vcf.gz

2 對所有樣品的gvcf文件進(jìn)行合并

用到的GATK模塊為CombineGVCFs, -Xmx400g -Xms400g 指定最大內(nèi)存，-V 每個樣品的gvcf文件 -O 指定輸出文件名稱

gatk --java-options "-Xmx400g -Xms400g -XX:+UseSerialGC"
CombineGVCFs -R reference.fasta
-V example1.g.vcf.gz -V example2.g.vcf.gz -V example3.g.vcf.gz
-O combine.g.vcf.gz

3 群體變異檢測

用到的GATK模塊為GenotypeGVCFs, -V gvcf文件 -O 指定輸出的vcf文件

gatk --java-options "-Xmx400g -Xms400g -XX:+UseSerialGC"
GenotypeGVCFs -R reference.fasta
-V combine.g.vcf.gz
-O combine.call.vcf.gz

4 提取SNP

用到的GATK模塊為SelectVariants, -select-type 指定提取類型

gatk --java-options "-Xmx400g -Xms400g -XX:+UseSerialGC"
SelectVariants -R reference.fasta
-V combine.call.vcf.gz
-select-type SNP
-O combine.SNP.vcf.gz

4 過濾SNP

用到的GATK模塊為VariantFiltration, --filter-expression 設(shè)置過濾參數(shù)及閾值

--filter-name 為每個位點(diǎn)添加標(biāo)記，大于過濾閾值的標(biāo)記為PASS，否則標(biāo)記為Filter

gatk --java-options "-Xmx400g -Xms400g -XX:+UseSerialGC"
VariantFiltration -R reference.fasta
-V combine.SNP.vcf.gz
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0
|| SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"
--filter-name "Filter"
-O combine.SNP.filter.vcf.gz

5 提取過濾好的SNP

用到的GATK模塊為SelectVariants，--exclude-filtered 去除被過濾掉的SNP，即標(biāo)記為Filter的位點(diǎn)

gatk --java-options "-Xmx400g -Xms400g -XX:+UseSerialGC"
SelectVariants -R reference.fasta
-V combine.SNP.filter.vcf.gz
--exclude-filtered
-O combine.SNP.filtered.vcf.gz

###############################################################
####### #######
####### 3 變異過濾 #######
####### 20220109 #######
###############################################################

VCF文件：example.vcf

1 去除Indel附近7bp范圍的SNP

-g 設(shè)置距離，單位bp， -O 指定輸出文件格式 v表示未壓縮的vcf文件

注：此步需要vcf文件中包含indel，過濾完后將indel去除

bcftools filter -g 7 -O v -o 1.vcf example.vcf

2 去除10bp范圍內(nèi)有大于3個SNP的SNP cluster

1 為vcf文件建立索引

gatk IndexFeatureFile -F 1.vcf

2 為符合要求的位點(diǎn)添加標(biāo)記SNP cluster，-cluster 指定SNP數(shù)目 -window 指定窗口范圍，單位bp

gatk VariantFiltration -V 1.vcf -cluster 3 -window 10 -O 2.vcf

3 去除SNP cluster

gatk --java-options "-XX:+UseSerialGC"
SelectVariants -V 2.vcf
-select "FILTER == SnpCluster"
--invertSelect
-O 2.Filtered.vcf

3 將質(zhì)量值小于20，覆蓋的reads數(shù)小于5的基因型設(shè)為miss即./.

vcftools --vcf 2.Filtered.vcf
--minDP 5 --minGQ 20
--recode --recode-INFO-all
--out 3.Filtered

4 去除非二等位基因、缺失率大于90%，次要等位基因頻率小于5%的位點(diǎn)

vcftools --vcf 3.Filtered.recode.vcf
--max-allele 2 --min-allele 2
--max-missing 0.9 --maf 0.05
--recode --recode-INFO-all
--out 4.Filtered

5 過濾連鎖不平衡位點(diǎn)

1 轉(zhuǎn)換vcf文件為ped格式

plink --vcf 4.Filtered.recode.vcf
--allow-extra-chr --recode
--out 5.ld

2 將map文件中的第二列改為每個snp位點(diǎn)特有的標(biāo)識符

awk '{print 1"_"3"\t"$4}' 5.ld.map > 5.ld1.map

3 修改相應(yīng)的ped文件的名稱

mv 5.ld.ped 5.ld1.ped

4 計(jì)算位點(diǎn)之間的ld值, --indep-pairwise 設(shè)置窗口、步長、r2值

plink --file 5.ld1
--allow-extra-chr
--indep-pairwise 50 10 0.2
--out ld

5 抽出低連鎖的位點(diǎn)，--extract 抽取滿足要求的位點(diǎn)，--make-bed 轉(zhuǎn)換為bed格式

plink --file 5.ld1
--allow-extra-chr
--extract ld.prune.in
--make-bed --out filter_ld

6 將bed文件轉(zhuǎn)為vcf文件

plink --bfile filter_ld
--allow-extra-chr
--recode vcf-iid --out filter_ld

###############################################################
####### #######
####### 4 系統(tǒng)樹、Admixture和PCA #######
####### 20220109 #######
###############################################################

VCF文件：example.vcf
BED文件：example.bed

系統(tǒng)樹

1 將SNP轉(zhuǎn)為fasta格式，用到的軟件為vcf2phylip.py(https://github.com/edgardomortiz/vcf2phylip)

python vcf2phylip.py -f -i example.vcf

2 megacc 構(gòu)建NJ樹

infer_NJ_nucleotide_NOld.mao為參數(shù)文件，常用參數(shù)如下：

[ AnalysisSettings ]

Analysis = Phylogeny Reconstruction

Scope = All Selected Taxa

Statistical Method = Neighbor-joining #指定建樹方法

Phylogeny Test = ====================

Test of Phylogeny = Bootstrap method

No. of Bootstrap Replications = 1000 #指定Bootstrap次數(shù)

Substitution Model = ====================

Substitutions Type = Nucleotide

Model/Method = p-distance #指定模型

Number of Threads = 4 #指定線程數(shù)目

-a 指定參數(shù)文件，-d指定序列文件

megacc -a infer_NJ_nucleotide_NOld.mao -d example.fasta -o example_out

3 IQTREE 構(gòu)建ML樹

-s 指定序列文件，-o 指定外類群名稱， -T 指定線程數(shù)目，-B 指定Bootstrap次數(shù)，-m MFP 自動搜索最適合模型并建樹

iqtree -s example.fasta -o acoe -T 10 -B 1000 -m MFP

Admixture

1 依次計(jì)算K=2到K=20

for k in {2..20} ;do admixture --cv=10 example.bed {k}.out; done

2 處理admixture結(jié)果

1 提取出Q文件中的個體順序

awk '{print $2}' example.fam > ID

2 給Q文件添加個體id

for i in *Q ;do paste ID i.id ;done

3 按照sortID文件中的指定順序，給Q文件排序

for i in *id ;do python3 sortQ.py sortID i.sort ;done

3 繪圖

a <- read.table("sort.5k",row.names = 1)
b <- t(as.matrix(a))
b.reverse <- b[,324:1]
p <- barplot(b.reverse,space = 0,border = NA,axisnames = F,
col = c("gold","firebrick","DarkGoldenrod4","IndianRed1",))
axis(side = 1, p, labels = F, at = c(0:324))
labs <- colnames(b.reverse)
text(cex=.3, x=p-0.25,y=-0.05,labs,xpd=TRUE,srt=45,adj = c(1,1))

PCA

--pca 324 計(jì)算324個PC的值，324為輸入文件中的個體數(shù)

plink --allow-extra-chr --bfile example --pca 324