浅谈如何在优麒麟22.04中使用Eigenstrat和Plink工具生成类23andMe格式原始数据

在使用此教程前，请先在Oracle VM VirtualBox或其他虚拟机中安装好优麒麟22.04版系统，或其对应源系统Ubuntu 22.04。

一、下载1240K数据包

https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-datahttps://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data

这里以基因数据包V50.0.p1的1240K数据为例。

二、下载并部署EIGENSTRAT

虽然GitHub上也有对应的非官方Windows版程序，但无论32位还是64位，都存在文件的输出大小超过2GB导致无法转化的问题，Windows子系统版的Linux（WSL）也有一定局限性，因此这一块需要在Linux系统本体或虚拟机中执行。

https://data.broadinstitute.org/alkesgroup/EIGENSOFT/EIG-6.1.4.tar.gzhttps://data.broadinstitute.org/alkesgroup/EIGENSOFT/EIG-6.1.4.tar.gz

或者使用wget+tar指令来下载安装。该文件里的程序包不用再编译安装（官方GitHub源下载的EIG需要手动编译安装）。由于优麒麟系统缺少一些依赖来运行convertf，因此还需要到Ubuntu官网下载：
https://packages.ubuntu.com/https://packages.ubuntu.com/

测试的过程中我安装了以下系统程序依赖后，convertf方可正常运行：

gcc-6-base_6.4.0-17ubuntu1_amd64.deb
gcc-8-base_8.4.0-1ubuntu1~18.04_amd64.deb
libc6_2.27-3ubuntu1.5_amd64.deb
libgcc1_8.4.0-1ubuntu1~18.04_amd64.deb
libgfortran3_6.4.0-17ubuntu1_amd64.deb
libquadmath0_8.4.0-1ubuntu1~18.04_amd64.deb

将convertf复制粘贴到数据包所在文件夹中，然后在对应路径新建空文本（重命名为：par.EIGENSTRAT.PACKEDPED），文档内容如下：

genotypename:    v50.0.p1_1240K_public.geno
snpname:         v50.0.p1_1240K_public.snp
indivname:       v50.0.p1_1240K_public.ind
outputformat:    PACKEDPED
genotypeoutname: v50.0.p1_1240K_public.bed
snpoutname:      v50.0.p1_1240K_public.bim
indivoutname:    v50.0.p1_1240K_public.fam

在要执行的路径下右击打开终端，输入

./convertf -p par.EIGENSTRAT.PACKEDPED

三、下载并部署Plink

PLINK 1.9https://www.cog-genomics.org/plink/1.9/解压plink在该路径后下右击打开终端，输入

./plink --bfile v50.0.p1_1240K_public --export vcf --out v50.0.p1_1240K_public

四、裁剪刚生成的VCF文件

打开“.ind”或“.fam”文件并对照信息的“行数+9”来切割需要的数据，使结果对应的列只有一列（因为转化原始数据时Plink不支持VCF文件种包含多列结果）。接着使用plink，也可以在文件生成后删除数据中前面所有以“#”开头的行。

比如，数据Denisova_published.DG在indiv文件中对应第501行，因此切出VCF文件前9列和第“501+9=510”列，并转换：

cut -f 1-9,510 v50.0.p1_1240K_public.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23 --out ./RawData/Denisova_published.DG.txt

或者去除所有包含“#”符号的行：

cut -f 1-9,510 v50.0.p1_1240K_public.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/Denisova_published.DG.txt

如果一次提取多个数据，那么分两次裁切VCF文件数据，示例如下：

cut -f 1-9,307-309,313-314,372 v50.0.p1_1240K_public.vcf > test6.vcfcut -f 1-9,10 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Sardinian-3.DG.txtcut -f 1-9,11 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Yoruba-3.DG.txtcut -f 1-9,12 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Karitiana-3.DG.txtcut -f 1-9,13 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Papuan-15.DG.txtcut -f 1-9,14 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Ju_hoan_North-4.DG.txtcut -f 1-9,15 test6.vcf > temp.vcf
./plink --vcf temp.vcf --recode 23
sed '/^#/d' plink.txt > ./RawData/B_Han-3.DG.txt

最后删除生成的中间文件。

注意：

1. 也可以使用 bcftools query -f '%ID\t%CHROM\t%POS[\t%TGT]\n' 指令来转化VCF文件，接着在用sed指令转化之后切裁。

2. Plink转化出的类23andMe格式原始数据的Y染色体和线粒体MT的检出位点只显示单个字母而非成对出现，可能会影响少数特殊的使用情况。

3. 在用Plink转23andMe格式前，如果VCF文件的样本名中包含下划线“_”，会导致转化出错，需先想办法转化或去除。

4. Reich数据包里有少数数据的样本名备注有误，需慎用，或者另寻科研论文对应的基因数据包。

5. 也可以在convertf转化步骤后先使用 plink --bfile v50.0_1240k_public --keep list_fam.txt --make-bed --out TEST_1240k 指令来实现提取所需要的样本，其中 list_fam.txt 为从“fam”后缀名文件复制出的所需样本所在行的集合；然后再依次转化为VCF与23andMe格式，以大幅节省计算机工作量。（此次新增内容）

扩展阅读：

https://zhuanlan.zhihu.com/p/73046966https://zhuanlan.zhihu.com/p/73046966https://github.com/DReichLab/EIGhttps://github.com/DReichLab/EIGhttps://reich.hms.harvard.edu/software/InputFileFormatshttps://reich.hms.harvard.edu/software/InputFileFormatshttps://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-datahttps://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data