前言
當歲月靜好，一切就等著到“曼哈頓”相見的時候
“BUG”就悄悄地來了！～
如暴風(fēng)雨前一樣安靜。
若數(shù)據(jù)安好，便是晴天～
當一切如前話獲得一個又一個的關(guān)聯(lián)結(jié)果后，卡殼的事情出現(xiàn)了！
當我使用keep保留希望留住想要的亞組的時候，居然報錯了?。ㄊ堑?-23號染色體，只有15號有問題）
提取完后居然沒有剩下的樣本??？
反正我也有22條染色體了～ 不如......
騷年！忘記騷想法！科學(xué)精神都去哪里了?。?！
于是就要開始查原因了。
人話：結(jié)果是簡單的，但是查找原因的過程是痛苦的～
正如：相愛是簡單的，但是相處的日子是痛苦的一樣。
第一步：一開始，我并不知道Individual ID是不是真的不存在（當然不存在的可能性不大～）
但是本著科學(xué)的精神，我們還是在需要提取的txt文件（--keep中間的那個文件）中挑一個individual ID出來，
grep -nw 142120126 chr19.ped?
IID說：你愛或不愛，我就在哪里～
那么，究竟是什么阻擋著我們的康莊大道？
第二步：為什么其他染色體能夠提取出來，就這個不行？
是否有天生異品之處？
我們先看看文件長什么樣子！
打開15號（不能提取的），19號和22號（能提取的）文件，看長什么樣子（文件較大，打開需要一點時間～ 注意加head?。。?！否則普通電腦handle不了?。?/p>
awk '{print $1,$2}' chr15.ped |head?
awk '{print $1,$2}' chr19.ped|head
awk '{print $1,$2}' chr22.ped |head
結(jié)果：
第三步：把所有不同的??給找出來！
你想想，就看前幾行就已經(jīng)有不一樣了，中間的呢？后面的呢？
相信我，你是不可能全部打開，然后一個一個對比的?。。?！
（有聽過愚公移山嗎？）
我們首先讀入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）認識文件名）
查看一下，我們關(guān)心的第二列（Individual ID，我們想保留的）是否不一樣
#首先配對
tmp=match(chr15fam[,2], chr20fam[,2])
#然后查看沒有配對的元素個數(shù)（我們的目的是查不同啊～）
length(which(is.na(tmp)))
結(jié)果：
[1] 2
沒想到還真有?。?5號染色體有兩個不同于其他染色體的存在?。?/p>
究竟是哪個?。?！
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，發(fā)現(xiàn)了某些individual ID缺失了，導(dǎo)致和FamilyID的配對亂了！
所以我們盡管看看是否如此：
查看15號染色體第274行（減一行）的樣子
chr15fam[274,]
結(jié)果：
? ? ? ? ? V1? ? ? ? V2 V3 V4 V5 V6
274 famid274 142120304? 0? 0? 0 -9
對比15號染色體第275行
chr20fam[275,]
好的，發(fā)現(xiàn)第二行的Individual ID是一樣的，但是Family ID確實不一樣！
? ? ? ? ? V1? ? ? ? V2 V3 V4 V5 V6
275 famid275 142120304? 0? 0? 0 -9
這個就是為什么之前的keep完全提取不了數(shù)據(jù)的原因?。?！
我們要知道，--keep中間的文件是有兩行的，是配對的！
這個染色體亂套了??！
是的，一子錯滿盤皆落錯！（這是原來數(shù)據(jù)集的問題啦～）
第四步：我們通過科學(xué)的方法查找了問題所在了，接下來就是解決問題了！～
#我們首先在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)
#可省略（不過習(xí)慣上都要看看數(shù)據(jù)讀入成不成功，不成功的話，后面的分析都是不靠譜的！～）
scfam[1:2,]
chr15[1:3,]
?dim(chr15)
dim(chr20)
control15[1:3,]
#我們再確認一下，第二行是否有差異，差異的是哪個（上述第三步重復(fù)）
tmp=match(chr15[,2], chr20[,2])
length(which(is.na(tmp)))
which(is.na(tmp))
chr20[275,]
chr20[1211,]
Terminal界面：
大招來了！～
設(shè)定一個變量（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號染色體需要提取的數(shù)據(jù)（scforchr15）：
對于15號染色體的行，根據(jù)scinchr15保留，對于15號染色體的列，全部第一，二列全部保留。
scforchr15= chr15[scinchr15, 1:2]
#scforchr15[1:3,]
#scinchr15[1:10]
#輸出數(shù)據(jù)，用于--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)
后記
告訴大家一個不幸的消息：盡管我經(jīng)過如此謹慎的計算后，該染色體并沒有我們想要的SNP
......
