使用SNP数据计算IBD、PCA、系统树鉴别亲缘关系

小故事

接到的小任务，对于大佬们来说可能不难，对于新手，却是要抓破脑袋。
两个养殖户的牛群跑到一个山头，牛群混了，有一头小牛不知道是谁家的，都坚称是自己的。他们自己掏钱找公司分别测序了两头母牛和小牛，数据送到导师手里。

数据预处理

1.得到的fastq文件，首先去除质量低的位点（QC）。
软件：NGSQCToolkit
2.序列比对，sentieon生成sort.bam，sambamba建立索引。
软件：sentieon、sambamba.
参考基因组需要另外下载。
3. 去冗余
软件：sambamba

sambamba markdup -r -t 10 --tmpdir ${prefix}/$index $prefix/$index/${index}_sort.bam $prefix/$index/${index}_mkdup.bam

–tmpdir 更换缓存目录，避免公共目录太满影响运行。
4.重点来了，合并三个mkdup.bam文件，这个时候合并bam文件，一直到后面生成vcf,软件都会识别成3个样本。当初合并vcf文件再执行其他操作时出现了不少问题。

//可以把mkdup.bam文件移动到同一目录下，用*号代表所有前缀。
samtools merge beefmaster_merge.bam *.mkdup.bam
sambamba index beefmaster_merge.bam

5.call SNP.
听讲座的时候老师一直说"sleep"，我一直没反应过来，原来说的是SNP
sentieon大法好，就是快。代码是老师的，不能随意公布。

6.过滤SNP.
用gatk过滤

for iterm in /public/home/********/dabieshan_merge
do
index=$(basename $iterm)
prefix=$(dirname $iterm)
cd $prefex
module load GATK
bsub -J gatk -e %J.gatk.err "
gatk SelectVariants \
    -select-type SNP \
    -V ${iterm}.snp.vcf.gz \
    -O ${iterm}.select_snp.vcf.gz
gatk VariantFiltration \
    -V ${iterm}.select_snp.vcf.gz \
   --filter-expression 'QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 4.0 || ReadPosRankSum < -8.0' --filter-name 'Filter' \
    -O $prefix/${index}/${index}_snps_filtered.vcf.gz"
done

我自己写的随意丢出来了。用的学校集群，有共享软件库，module load 引入软件。

计算IBD

软件：Plink v-1.9

plink --vcf dabieshan_merge_snps_filtered.vcf.gz --recode --out dabieshan_merge_snps_filtered --const-fid --allow-extra-chr
// --const-fid 添加类群名，转换格式
plink --allow-extra-chr --file dabieshan_merge_snps_filtered --noweb --make-bed --out dabieshan_merge_snps_filtered
// 得二进制文件
plink --allow-extra-chr --file dabieshan_merge_snps_filtered --noweb --maf 0.15 --geno 0.03 --indep-pairwise 1000 50 0.2 --genome --out dabieshan_merge_snps_filtered
// 过滤
// 其实可以也可以不用生成这么多文件，直接计算IBD，这么多文件只是方便后面计算PCA和建树
plink --vcf dabieshan_merge_snps_filtered.vcf.gz  --const-fid --allow-extra-chr --maf 0.15 --geno 0.03 --indep-pairwise 1000 50 0.2 --genome --out IBD_result

调用–genome命令就是生成。genome文件，里面就是IBD结果

           FID1           IID1           FID2           IID2 RT    EZ      Z0      Z1      Z2  PI_HAT PHE       DST     PPC   RATIO
  FDSW202195227  FDSW202195227  FDSW202195228  FDSW202195228 UN    NA  1.0000  0.0000  0.0000  0.0000  -1  0.594297  0.0000  1.0512
  FDSW202195227  FDSW202195227  FDSW202195229  FDSW202195229 UN    NA  0.7350  0.2650  0.0000  0.1325  -1  0.638245  0.0000  1.4630
  FDSW202195228  FDSW202195228  FDSW202195229  FDSW202195229 UN    NA  0.3546  0.6454  0.0000  0.3227  -1  0.716348  0.0003  1.8117

PI_HAT值越大，亲缘关系越近，PI_HAT=1/2*Z1 + Z2，为1则可能为同卵双胞胎或者重复测序样本。Z0表示样本IID1与IID2的0个染色体组含有相同变异的频率,同理，z1/z2表示两个样本1个/2个染色体组含有相同变异频率。
这个IBD结果表明，227号与228号关系远，为两头母牛的测序号。
228与229关系更近，两样本Z1为0.6454，初步断定228为小牛229的母亲。当然还需PCA、系统树的辅助验证。
**实践证明：合并vcf文件计算得到的PI_HAT大的离谱，合并bam文件更合理。

计算pca

上一步以及生成二进制文件，用bfile指定文件。

plink --allow-extra-chr --bfile dabieshan_merge_snps_filtered --noweb --maf 0.15 --geno 0.03 --indep-pairwise 1000 50 0.2 --pca 5 --out dabieshan_merge_snps_filtered

–maf 0.15 --geno 0.03 --indep-pairwise 1000 50 0.2为过滤条件
–pca 会生成pca结果文件.eigenval和.eigenvec。数字5表示需要的主成分数量，默认为20，我这里只有3个样本，不能找得出这么多成分。文件eigenval中记录了主成分的特征值，特征值得和约等于样本数，可根据这个设置pca后的数值。
再用.eigenvec文件作图。在文件头行插入一列
原文件

FDSW202195227 FDSW202195227 0.803103 -0.147283 -0.57735
FDSW202195228 FDSW202195228 -0.529102 -0.621866 -0.57735
FDSW202195229 FDSW202195229 -0.274001 0.769149 -0.57735

插入头行

FID IID pc1 pc2 pc3
FDSW202195227 FDSW202195227 0.803103 -0.147283 -0.57735
FDSW202195228 FDSW202195228 -0.529102 -0.621866 -0.57735
FDSW202195229 FDSW202195229 -0.274001 0.769149 -0.57735

再用R作图

// 画图，样式大小和颜色可调整。pc1的比例大，则坐标刻度越大。
plot(dabieshan$pc1,dabieshan$pc2, pch=c(1,2,3), col=c(1,2,3) , cex=3, main="PCA",xlab="pc1 (46%)",ylab="pc2 (28%)")
// 图注，可调整图注位置
legend("topright",c("FDSW202195227","FDSW202195228","FDSW202195229"), pch=c(16,17,18), col=c(1,2,3))

距离越近，则越聚在一起。任然是228和229关系近。

系统发育树

Github搜索vcf2phylip.py，一个开源脚本，持续更新中，且有详细用法介绍。把vcf文件转成phy文件，一个MEGAX可识别的格式。复制代码后改下权限即可用。
导入MEGAX，转成meg格式，即可作NJ树图，由于样本较少，选择参数时方法选为none，否则会报错。

任然是228和229关系更近。

计算SNP分布情况。

每个样本单独得到vcf文件。根据SNP位点及基因型，作Venn图。导师自己写R脚本做的，不便分享。

后记

导师说结果已经发给两位养殖户，小牛是否寻得亲妈就不得而知了。测序近4000，我们的分析费用为0，不得不说当代农民越来越有钱了。能参与这么有意思的一件事还是挺开心的。