4.多个同样行列式文件的合并
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
这是一个样本的基因表达矩阵文件，但我们在处理实际数据的时候，往往不止有一个样本，所以对多个文件进行合并，得到一个包含全部样本基因表达情况的表达矩阵。
原始数据：
EnsEMBL_Gene_ID sample1
ENSG00000000003 2.2938981848
ENSG00000001167 61.6255191139
ENSG00000005471 0.5565311150
ENSG00000066629 32.5281120565
ENSG00000154258 0.1778082057
ENSG00000154262 0.1057582795
ENSG00000154263 0.2211719508
ENSG00000154265 4.3831737628
ENSG00000154269 0.0000000000
目标数据：
| EnsEMBL_Gene_ID | sample1 | sample2 | sample3 |
|---|---|---|---|
| ENSG00000000003 | 2.2938981848 | 0.1778082057 | 2.2211719508 |
| ENSG00000005471 | 1.1778082057 | 0.1057582795 | 4.3831737628 |
| ENSG00000154258 | 4.3831737628 | 0.2211719508 | 2.2938981848 |
| ENSG00000154269 | 2.2938981848 | 4.3831737628 | 2.2938981848 |
## 首先使用测试数据进行处理
### 创建测试代码
#新建tmp.sh文件，并在其中输入：
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
#使用perl命令调用
perl -e 'system("bash tmp.sh > $_.txt") foreach a..z'
生成了a-6共26个的counts文件，每个文件都包含两列，第一列为基因，第二列为表达量
使用R语言进行处理
#使用awk对全部文件进行合并
awk '{print FILENAME"\t"$0}' * | grep -v gene_ID > tmp.txt
# 使用R语言的reshape包进行处理
## 读入文件
a = read.table('tmp.txt',sep='\t', stringsAsFactors = F)
colnames(a) = c('Samples','Gene_ID','Count')
## 格式转换
fpkm <- dcast(a, formula = Gene_ID~Samples)
==dcast函数介绍==
dcast 将长数据转换成想要的任何宽型数据，并以数据框形式返回 dcast(data, formula, fun.aggregate = NULL, ..., margins = NULL, subset = NULL, fill = NULL, drop = TRUE, value.var = guess_value(data)) data 需要转换的数据，格式为数据框data.frame formula 用于转换的公式，格式为行：行变量～列变量，在这里，我们的目标数据格式应该是Gene_ID为行变量，而样本名为列变量
使用python进行转换
python中也有一个从长数据转换到宽数据的方法
长型数据：变量中至少有一个变量中的元素存在严重重复循环的情况 宽型数据：各变量的值不存在重复循环的情况也无法归类
import pandas as pd
#读入文件，添加参数header=None,防止误将第一行读为行名
df = pd.read_table('tmp.txt',header=None)
#给数据框添加列名
df.columns = ['Samples','Gene_ID','Count']
#使用数据透视进行转换，index为行变量，columns表示列变量，values为值
rpkm<-df.pivot_table(index=["Gene_ID"],
columns=["Samples"],
values=["Count"])
从ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tar下载数据，共56个文件
首先将tar文件解压
tar -xf GSE48213_RAW.tar
将全部样本的表达值文件合并为一个
ls *.txt.gz | cut -d'.' -f 1| sort -u | while read id; do zcat ${id}.txt.gz | awk '{print "'${id}'\t"$0}' | grep -v EnsEMBL_Gene_ID >> GSM.txt; done
使用R语言的reshape2的dcast函数进行转换
a = read.table('GSM.txt',sep='\t',stringsAsFactors=F,fill=T)
fpkm <- dcast(a,formula=Gene_ID~Samples)