前言
當(dāng)歲月靜好，一切就等著到“曼哈頓”相見的時(shí)候
“BUG”就悄悄地來了吞加！～
如暴風(fēng)雨前一樣安靜裙犹。
若數(shù)據(jù)安好，便是晴天～
當(dāng)一切如前話獲得一個(gè)又一個(gè)的關(guān)聯(lián)結(jié)果后衔憨，卡殼的事情出現(xiàn)了叶圃！
當(dāng)我使用keep保留希望留住想要的亞組的時(shí)候，居然報(bào)錯(cuò)了＜肌（是的1-23號(hào)染色體掺冠，只有15號(hào)有問題）
提取完后居然沒有剩下的樣本Ｂ氲场德崭？
反正我也有22條染色體了～ 不如......
騷年！忘記騷想法揖盘！科學(xué)精神都去哪里了Ｃ汲！兽狭！
于是就要開始查原因了憾股。
人話：結(jié)果是簡單的，但是查找原因的過程是痛苦的～
正如：相愛是簡單的箕慧，但是相處的日子是痛苦的一樣服球。
第一步：一開始，我并不知道Individual ID是不是真的不存在（當(dāng)然不存在的可能性不大～）
但是本著科學(xué)的精神颠焦，我們還是在需要提取的txt文件（--keep中間的那個(gè)文件）中挑一個(gè)individual ID出來斩熊，
grep -nw 142120126 chr19.ped?
IID說：你愛或不愛，我就在哪里～
那么伐庭，究竟是什么阻擋著我們的康莊大道粉渠？
第二步：為什么其他染色體能夠提取出來分冈，就這個(gè)不行？
是否有天生異品之處渣叛？
我們先看看文件長什么樣子丈秩！
打開15號(hào)（不能提取的），19號(hào)和22號(hào)（能提取的）文件淳衙，看長什么樣子（文件較大蘑秽，打開需要一點(diǎn)時(shí)間～ 注意加head！ｓ锱省３ι！否則普通電腦handle不了Ｑヵ恕）
awk '{print $1,$2}' chr15.ped |head?
awk '{print $1,$2}' chr19.ped|head
awk '{print $1,$2}' chr22.ped |head
結(jié)果：
第三步：把所有不同的??給找出來！
你想想绝葡，就看前幾行就已經(jīng)有不一樣了深碱，中間的呢？后面的呢藏畅？
相信我敷硅，你是不可能全部打開，然后一個(gè)一個(gè)對(duì)比的Ｓ溲帧＝时摹！榜旦！
（有聽過愚公移山嗎幽七？）
我們首先讀入15號(hào)染色體和20號(hào)染色體（隨便選的，一個(gè)能跑的就可以了）
打開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）認(rèn)識(shí)文件名）
查看一下，我們關(guān)心的第二列（Individual ID亭饵，我們想保留的）是否不一樣
#首先配對(duì)
tmp=match(chr15fam[,2], chr20fam[,2])
#然后查看沒有配對(duì)的元素個(gè)數(shù)（我們的目的是查不同啊～）
length(which(is.na(tmp)))
結(jié)果：
[1] 2
沒想到還真有Ｐ菖肌（15號(hào)染色體有兩個(gè)不同于其他染色體的存在！）
究竟是哪個(gè)９佳颉Ｌざ怠词顾！
which(is.na(tmp))
[1]? 275 1211
我們來看一個(gè)全相！
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的配對(duì)亂了！
所以我們盡管看看是否如此：
查看15號(hào)染色體第274行（減一行）的樣子
chr15fam[274,]
結(jié)果：
? ? ? ? ? V1? ? ? ? V2 V3 V4 V5 V6
274 famid274 142120304? 0? 0? 0 -9
對(duì)比15號(hào)染色體第275行
chr20fam[275,]
好的疹尾，發(fā)現(xiàn)第二行的Individual ID是一樣的上忍，但是Family ID確實(shí)不一樣！
? ? ? ? ? V1? ? ? ? V2 V3 V4 V5 V6
275 famid275 142120304? 0? 0? 0 -9
這個(gè)就是為什么之前的keep完全提取不了數(shù)據(jù)的原因Ｄ杀尽Ｇ侠丁！
我們要知道繁成，--keep中間的文件是有兩行的吓笙，是配對(duì)的！
這個(gè)染色體亂套了＝硗蟆面睛！
是的，一子錯(cuò)滿盤皆落錯(cuò)Ｗ鸢帷（這是原來數(shù)據(jù)集的問題啦～）
第四步：我們通過科學(xué)的方法查找了問題所在了叁鉴，接下來就是解決問題了！～
#我們首先在R中讀入原來用在--keep中的txt文件
scfam=read.table("SCfam.txt",header=F,stringsAsFactors=F)
#在terminal中提取15號(hào)染色體ped文件前兩行（Family ID和Individual ID）毁嗦，保存為一個(gè)txt文件亲茅。
awk '{print $1,$2}' chr15.ped > ~seedson/Desktop/file/Chr15IDs.txt
awk '{print $1,$2}' chr20.ped > ~seedson/Desktop/file/Chr20IDs.txt （用于后面的確認(rèn)而已）
在R中讀入15號(hào)染色體的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,]
#我們?cè)俅_認(rèn)一下，第二行是否有差異腔长，差異的是哪個(gè)（上述第三步重復(fù)）
tmp=match(chr15[,2], chr20[,2])
length(which(is.na(tmp)))
which(is.na(tmp))
chr20[275,]
chr20[1211,]
Terminal界面：
大招來了捞附！～
設(shè)定一個(gè)變量（scinchr15）巾乳，把需要提取的Individual ID和全部的15號(hào)染色體的Individual ID配對(duì)起來（第二列）
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]
對(duì)于15號(hào)染色體需要提取的數(shù)據(jù)（scforchr15）：
對(duì)于15號(hào)染色體的行，根據(jù)scinchr15保留鸟召，對(duì)于15號(hào)染色體的列胆绊，全部第一，二列全部保留欧募。
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，腰不疼啦,腿不痛啦..上六樓啊也有勁勒
對(duì)照組也是同樣的處理（其他亞組也是一樣了）
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)
后記
告訴大家一個(gè)不幸的消息：盡管我經(jīng)過如此謹(jǐn)慎的計(jì)算后种冬，該染色體并沒有我們想要的SNP
......
