前言
当岁月静好，一切就等着到“曼哈顿”相见的时候
“BUG”就悄悄地来了！～
如暴风雨前一样安静。
若数据安好，便是晴天～
当一切如前话获得一个又一个的关联结果后，卡壳的事情出现了！
当我使用keep保留希望留住想要的亚组的时候，居然报错了！（是的1-23号染色体，只有15号有问题）
小白幼小而又脆弱的心灵！
提取完后居然没有剩下的样本！？
反正我也有22条染色体了～ 不如......
骚年！忘记骚想法！科学精神都去哪里了！！！
于是就要开始查原因了。
人话：结果是简单的，但是查找原因的过程是痛苦的～
正如：相爱是简单的，但是相处的日子是痛苦的一样。
第一步：一开始，我并不知道Individual ID是不是真的不存在（当然不存在的可能性不大～）
但是本着科学的精神，我们还是在需要提取的txt文件（--keep中间的那个文件）中挑一个individual ID出来，
grep -nw 142120126 chr19.ped
是存在的～
反复查找是否存在
IID说：你爱或不爱，我就在哪里～
那么，究竟是什么阻挡着我们的康庄大道？
第二步：为什么其他染色体能够提取出来，就这个不行？
是否有天生异品之处？
我们先看看文件长什么样子！
打开15号（不能提取的），19号和22号（能提取的）文件，看长什么样子（文件较大，打开需要一点时间～ 注意加head！！！！否则普通电脑handle不了！）
awk '{print $1,$2}' chr15.ped |head
awk '{print $1,$2}' chr19.ped|head
awk '{print $1,$2}' chr22.ped |head
结果：
看到了吗？15号染色体和其他的染色体不！一！样！
第三步：把所有不同的™️给找出来！
你想想，就看前几行就已经有不一样了，中间的呢？后面的呢？
相信我，你是不可能全部打开，然后一个一个对比的！！！！
（有听过愚公移山吗？）
我们首先读入15号染色体和20号染色体（随便选的，一个能跑的就可以了）
打开R（Rstudio也好，Terminal的R也罢）
读入chr15.fam文件
chr15fam=read.table("chr15.fam",header=F,stringsAsFactors=F)
读入chr20.fam文件
chr20fam=read.table("chr20.fam",header=F,stringsAsFactors=F)
我说过了，熟悉这些文件是什么，是很重要的！（GWAS分析-说人话（2）认识文件名）
查看一下，我们关心的第二列（Individual ID，我们想保留的）是否不一样
#首先配对
tmp=match(chr15fam[,2], chr20fam[,2])
#然后查看没有配对的元素个数（我们的目的是查不同啊～）
length(which(is.na(tmp)))
结果：
[1] 2
没想到还真有！（15号染色体有两个不同于其他染色体的存在！）
究竟是哪个！！！
which(is.na(tmp))
[1] 275 1211
我们来看一个全相！
chr20fam[notfound,]
V1 V2 V3 V4 V5 V6
275 famid275 142120304 0 0 0 -9
1211 famid1211 142120782 0 0 0 -9
我们从一开始的head，发现了某些individual ID缺失了，导致和FamilyID的配对乱了！
所以我们尽管看看是否如此：
查看15号染色体第274行（减一行）的样子
chr15fam[274,]
结果：
V1 V2 V3 V4 V5 V6
274 famid274 142120304 0 0 0 -9
对比15号染色体第275行
chr20fam[275,]
好的，发现第二行的Individual ID是一样的，但是Family ID确实不一样！
V1 V2 V3 V4 V5 V6
275 famid275 142120304 0 0 0 -9
这个就是为什么之前的keep完全提取不了数据的原因！！！
我们要知道，--keep中间的文件是有两行的，是配对的！
这个染色体乱套了！！
是的，一子错满盘皆落错！（这是原来数据集的问题啦～）
第四步：我们通过科学的方法查找了问题所在了，接下来就是解决问题了！～
#我们首先在R中读入原来用在--keep中的txt文件
scfam=read.table("SCfam.txt",header=F,stringsAsFactors=F)
#在terminal中提取15号染色体ped文件前两行（Family ID和Individual ID），保存为一个txt文件。
awk '{print $1,$2}' chr15.ped > ~seedson/Desktop/file/Chr15IDs.txt
awk '{print $1,$2}' chr20.ped > ~seedson/Desktop/file/Chr20IDs.txt （用于后面的确认而已）
在R中读入15号染色体的txt文件
chr15=read.table("Chr15IDs.txt",header=F,stringsAsFactors=F)
#可省略（不过习惯上都要看看数据读入成不成功，不成功的话，后面的分析都是不靠谱的！～）
scfam[1:2,]
chr15[1:3,]
dim(chr15)
dim(chr20)
control15[1:3,]
#我们再确认一下，第二行是否有差异，差异的是哪个（上述第三步重复）
tmp=match(chr15[,2], chr20[,2])
length(which(is.na(tmp)))
which(is.na(tmp))
chr20[275,]
chr20[1211,]
Terminal界面：
相信我，查BUG就是要反复确认，因为没有人会告诉你答案
大招来了！～
设定一个变量（scinchr15），把需要提取的Individual ID和全部的15号染色体的Individual ID配对起来（第二列）
scinchr15= match(scfam[,2], chr15[,2])
#length(which(is.na(scinchr15)))
#scinchr15[1:10]
#scfaminchr15=match(scfam[,1], chr15[,1])
#length(which(is.na(scfaminchr15)))
#scfaminchr15[1:10]
对于15号染色体需要提取的数据（scforchr15）：
对于15号染色体的行，根据scinchr15保留，对于15号染色体的列，全部第一，二列全部保留。
scforchr15= chr15[scinchr15, 1:2]
#scforchr15[1:3,]
#scinchr15[1:10]
#输出数据，用于--keep中间。
write.table(scforchr15, file="SCfam_forChr15.txt",quote=F,row.names=F,col.names=F)
看见右下角的Done，腰不疼啦,腿不痛啦..上六楼啊也有劲勒
配对成功
对照组也是同样的处理（其他亚组也是一样了）
chr15=read.table("Chr15IDs.txt",header=F,stringsAsFactors=F)
control15=read.table("controlfam.txt",header=F,stringsAsFactors=F)
controlscinchr15= match(control15[,2], chr15[,2])
scforchr15control= chr15[controlscinchr15, 1:2]
write.table(scforchr15control,file="control_forChr15.txt",quote=F,row.names=F,col.names=F)
后记
告诉大家一个不幸的消息：尽管我经过如此谨慎的计算后，该染色体并没有我们想要的SNP
......
原谅我的浮躁～
网友评论