前言

书接上一回：

keep并不是像如下的那么简答的！

不要被表象骗了

这个指令的关键在于中间的sampleID文件的准备

--keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column

这是官网上的说明，说明这个文件必须是要有两行的！！！第一行family IDs，第2行within-family IDs

一切都是为了准备这一个文件！

我们回顾一下fam文件：第一列和第二列还是Family ID、Individual ID

我们上一话，提取的是Individual ID，所以，现在就是要根据需要的Individual ID（上一话已提取好了），在fam文件中，提取染色体2fam文件（chr2.fam，说人话（3）已经准备好了）中需要的信息。

直接放大招：

#打开R studio

#寻找文件所在的文件夹位置

#本例子在这个目录（~/Documents/LungCancerSubtypes/）

#把之前提取的individual ID文件（SCLC_GENEVA_ID）拉进这个目录，不然后面操作不起来！～

#设置工作目录

setwd("~/Documents/LungCancerSubtypes/")

#fam也不是什么特别货色了，本质上也就是一个表格！

#用read.table 函数处理

但这里先做一个错误示范：

# fam = read.table ("chr2.fam", header = F, stringsAsFactors = F, sep = "\t")

#这里插入一个查错的方法，dim（非必要）

#dim(fam)

dim查错

#其实也就看看，数据的结构张什么样子的，多少行，多少列，心里有个B数

#一看这个结果只有一列就知道不对了

#因为熟悉fam文件的你，就知道这玩意不可能只有1列（不知道fam文件的，查看：https://www.jianshu.com/p/8ced7531b728），我尽力了。

正确的输入如下：

> fam = read.table ("chr2.fam", header = F, stringsAsFactors = F)

> dim(fam)

可见输入正确，这文件是用空格分开数值的：

正确的输入

#然后读入需要配对的ID文件（一开始拉到相同文件夹的玩意，上一回准备的文件）：

>SCLCID = read.table ("SCLC_GENEVA_ID.txt", header = F, stringsAsFactors = F)

> dim(SCLCID)

大招来了哦！

#定义匹配！！！match函数，少数向多数配对/靠拢/匹配

#我们需要的是根据SCLCID文件中的第一列（少数数据），提取fam文件（大文件）中的第二列

#因为第2列有齐SCLCID的数据（根据SCLCID的值，提取fam文件想要的数据）

temp1 = match(SCLCID[,1],fam[,2])

然后length函数一下，看看有多少数据（排查错误）

> length(temp1)

#插入知识点：

dim结果

#这里如果用dim，为什么会报错？

#这个temp1，只是一个向量（计算机编程定义：存放数据的地方，一维数组），不是数据集！dim没用的～（NULL）

#不过看一下这个向量，就知道配对所在的位置了哦！（这个会有其它方面的应用的）

> temp1[1:10]

#所以，其实就提取好了！

#我们要配对的行（temp1），该文件的第1，2列

>famSCLC=fam[temp1, 1:2]

#检查：

> dim(famSCLC)

#因为NA是一个不太受欢迎的存在，很多时候程序跑不出来识别不了。

> length(which(is.na(temp1)))

查看有没有NA

#可见，这个数据集还是有1个老鼠屎的

#接下来就是要出去了：没错，就是which配合叹号（逻辑运算： &，|，!。(与，或，非。)，这个复习，够全面了吧？）

> tempok=temp1[which(!is.na(temp1))]

#再次查看数据对不对（科研工作还是必须严谨的）

> length(tempok)

> tempok[75:80]

#再次提取，这才是我们想要的数据！

> famSCLC=fam[tempok, 1:2]

> dim(famSCLC)

#看最后5行张咋样

> famSCLC[75:80,]

#大功告成，写出文件（本例写到桌面上，文件命名（随意弄，就加上今天的日子吧），不要行名，也不要列名，空格作为分割符）

~~> write.table (famSCLC, file ="~/Desktop/famSCLC2019.txt", sep ="", row.names =FALSE, col.names =FALSE)~~ （下面有说明）

好了！我们写了足足两个文章，就是为了这个文件了！

famSCLC2019.txt

其实上面是给自己挖了个坑，因为没有用空格符"\t"分开两行，所以plink识别不了

解决方法一：

> write.table (famSCLC, file ="~/Desktop/famSCLC2019.txt", sep ="", row.names =FALSE, col.names =FALSE)

解决方法二：

重新读读入数据后再输出分开的文件，重命名或覆盖：

fam=read.table("famSCLC2019.txt",header=F, stringsAsFactors=F)

fam[1:2,]

write.table(fam, file="SCfam.txt",quote=F, row.name=F, col.names=F, sep="\t")

fam[1:2,]

请看下一文，正式使用keep指令吧～

GWAS分析-说人话（6）keep提取我们需要的表型数据 II