【PCA】scaffold等级的群体重测序可能遇到的坑
这阵子在做红树的群体重测序分析，一开始在PCA这个坎上遇到了一个网上普遍的教程都没有怎么提及的坑，加之下午有人发邮件给我，问了类似的问题，所以在这里写下这篇可能对大部分只有scaffold的参考基因在做pca分析时有帮助的文章。
在一开始，确定软件使用的时候，我在简书里找到西澳老哥lakeseafly的这一篇介绍gwas常用做pca的软件使用的文章，里面提到了三个软件，分别是gcta，eigensoft以及r里面一个叫vcfR的包。然后我就试了一下，结果发现gcta只设计了针对人的基因组的pca分析（即基因组组装到染色体等级并且少于22条染色体的才能使用），而vcfR里面有一个依赖包poppr这个包我安装不上，故而只剩下eigensoft这个软件可以用。
下载，安装之后，我就看了一下readme，然后被绕的有点晕，所以惯性的去请教鲍大哥怎么用。鲍大哥跟我讲了输入文件的格式以及需要注意的一些坑，然后我就按照他说的，开始了一波瞎折腾，然后就遇到了第一个坑，LD过滤报错？？？
具体的报错就是使用plink的-ld-pairwise这个参数时，它提示too much variants were pass，这个就有点蓝瘦了，然后我就用信鸽搜索，将报错打进去，大概知道了是由于snp莫得id导致，故而我再次用信鸽搜了一下怎么加id，然后就搜到了大神严涛的博客，在他那篇用emmax做gwas分析的博文中，我看到了黎明的曙光，原来gatk生成的vcf文件中是没有snp_id的，需要自己去加，具体怎么加，照搬大神的代码
plink --vcf test.vcf --maf 0.05 --geno 0.1 --recode vcf-iid --out test.maf0.05.int0.9 --allow-extra-chr
其中maf跟geno是常用的过滤条件，vcf-iid这个参数就是加入snp_id，由于没有达到染色体级别，所以需要用到allow-extra-chr这个参数
而后的话，过滤就正常的进行了下去，而之所以做这一步的意义在于过滤掉这些冗余的信息，加快计算的速度，而过滤的具体参数需要根据自己的物种特点（自交or杂交）去查类似的物种别人怎么过滤。
质控就算正常进行下去了，然后就是eigensoft的输入文件的准备了，在自己折腾了半天还是没搞对格式之后，我在qq上问了三楼的王师姐，师姐因为有事，所以让我去服务器上看一下她的输入文件，然后我就把师姐的输入文件down了一份下来，仔细观摩发现，我傻逼了，ped文件上次不知道哪个步骤被替换成另一种格式，但仍然以ped结尾，难怪搞半天没反应。意识到自己的傻逼了之后，我重新准备了输入文件，然后按照流程，迅速的得到了下一个坑。
只看到那smartpca检测完我的输入文件无误之后，开始了分析，而后，warning，bad chromosome，然后就一堆东西类似的出来之后，退出了。
看来在成为一名称职的帅哥前，还是需要接受很多磨难，然后我再次把报错内容交给了信鸽，然后就搜出了罪魁祸首，map文件中的染色体为0，嘤嘤嘤(╥╯^╰╥)，没有染色体水平的基因组就是个离了婚的臭男人。然后在众多搜索结果中，我发现了一位仁兄分析了他的报错经历。照着他说的，我把map文件第一列改为1，然后就成功运行了。
然后在那个外面天气我也不知道怎么样的下午，我又找到了宇宙实验媛的知乎，在她的知乎上，我看到了另一种输入文件格式。
在本人的不定期吐槽小文最后，稍微正经一下，把上面这些东西规整一下，别看得云里雾里。
第一步，使用plink进行过滤，并且加上snp的id，加法见上面的那道疑似的命令行，当然，更希望你们去看一下严涛大神的原文，他写得更好。
第二步，使用plink对加上id的vcf文件进行ld的过滤，得到prune.in跟out两个id文件，而后再使用plink的--extract选择in那个文件对加了id的vcf文件进行过滤。
第三步，将bed文件转为vcf文件，反正后面做其他分析可能还需要以vcf为输入文件，所以建议输出
第四步，eigensoft的输入文件准备，看下面：
输入文件有很多种形式可以选择，其中推荐两种，一种是ped文件作为输入格式，另一种是bed文件作为输入格式
首先是以ped文件为输入文件的格式
需要的三个文件分别是ped，map，pedind，需要注意的是这三种格式是怎么来的。
ped跟map直接由plink生成，pedind需要自己生成。
pedind其实就是ped文件的前六列，这个可以通过脚本轻易的实现。
将map文件第一列的0全部改为同一个大于0数字。
smartpca.perl -i myfile.ped -a myfile.map-b myfile.ind -o myfile.PCA -p myfile.plot -e myfile.evel -l myfile.log
这样就ok了。
第二种的输入文件是bed，bim，ind。
Bed，bim都是plink转换生成的，其中bim的第一列是染色体也一样需要改变一下
Ind由三列构成第一列为样本ID，第二列为性别（M，F，U），第三列为所属子群。
然后我们可以把输入文件写成一个文件，例如这样：
numoutevec是保留多少pca主成分
然后搞定之后就可以运行了，smartpca –p runningpca.conf > smartpca.log && echo"----smartpca----done"
这样就完成了pca的运算了。
然后，这里补充一下后面鲍大哥看了我的步骤提出来的几点修改建议，map文件的chr位置均为0是由于vcftools转换的时候不能识别，解决的方法可以是在转换的后面加上--chrom-map，然后给个对应的列表就加上去了，而这时添加snp_ID就并不是必要的。修改染色体名字就相当于强行把不同的scaffold变成一条scaffold。这样在计算的时候是有问题的。（由于我用别的软件算完pca结果其实跟修改map的结果差不多，所以这里其实是可以修改染色体的，但后面其他的运算可能就必须注意，所以大家得注意一下）
最后，如果不会脚本，也可以用plink进行pca运算
Plink –bfile myfile –pca 10 –out myfile.pca，简单，快捷，而且结果跟eigensoft差不多。
而后在看到牛奶的这篇知乎，觉得这个东西很有用，所以分享一下
当跑完pca之后，一般需要进行群体分层，分层的依据可以根据structure，pca以及nj_tree的综合判断，而至于pca要怎么判断，eigensoft里面其实内置了一个计算pca结果是否显著的工具，我们只要调用就行
twstats -t twtable -i myfile.pca.eval -o myfile.eigenvaltw.out
最后less看一下就行，p值小于0.05的都是显著
最后的最后，推广一下我们专栏跟公众号，知乎专栏【高通量测序技术】，我的知乎【生信小撰】，我们的微信公众号【高通量测序技术】，欢迎大家关注，也欢迎各位新手跟大神来互相交流。
如果大家遇到什么疑问，也欢迎将文件发送到这个邮箱：825526231@qq.com，只要是我知道的，我会尽力解答~
