其他
生物信息学技能面试题(第4题)-多个同样的行列式文件合并起来
相信用过htseq-count的朋友都知道,它是分开对每个样本计算所有的基因表达量,所以会生成一个个独立的文件,我用perl脚本模仿它的结果如下:
$ head a.txt
gene_1 178
gene_2 692
gene_3 486
gene_4 666
gene_5 395
gene_6 48
gene_7 926
gene_8 733
gene_9 660
gene_10 578
第一列是基因,第二列是该基因的counts值,共有a~z这26个样本的counts文件,需要合并成一个大的行列式,这样才能导入到R里面做差异分析,如果手工用excel表格做,当然是可以的,但是太麻烦,如果有500个样本,正常人都不会去手工做了,需要编程。
生成测试文件的代码如下:
#首先新建文件tmp.sh 输入这个代码:
perl
-
le '
{
print
"gene_$_\t"
.int
(
rand
(
1000
)
)
foreach
1.
.
99
}
'
## 然后用perl脚本调用这个tmp.sh文件:
perl
-
e 'system
(
" bash tmp.sh >$_.txt"
)
foreach a..z'
##这样就生成了a~z这26个样本的counts文件
用shell或者perl或者python,设置R语言都可以做,但是各有优缺点,而且如果每个样本的基因顺序并不一致,这时候你应该怎么做呢?
实际需求如下: 里面有56个文件(),需要合并成一个表达矩阵,来根据cell-line的不同,分组做差异分析。
paper是:
输出的表达矩阵,如下所示:
请务必做出此题,很重要!
先给一下shell结合R语言的做法:
## 首先在GSE48213_RAW目录里面生成tmp.txt文件
awk
'{print FILENAME"\t"$0}'
*
|grep
-
v EnsEMBL_Gene_ID >tmp.txt
## 然后把tmp.txt导入R语言里面用reshape2处理即可!
setwd(
'tmp/GSE48213_RAW/'
)
a
=
read.table(
'tmp.txt'
,sep
=
'\t'
,stringsAsFactors
=
F)
library(reshape2)
fpkm <
-
dcast(a,formula
=
V2~V1)