生物信息学技能面试题(第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语言都可以做,但是各有优缺点,而且如果每个样本的基因顺序并不一致,这时候你应该怎么做呢? 实际需求如下:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213 里面有56个文件(ftp://ftp.ncbi.nlm.nih.gov/geo/s ... pl/GSE48213_RAW.tar),需要合并成一个表达矩阵,来根据cell-line的不同,分组做差异分析。 paper是:https://www.ncbi.nlm.nih.gov/pubmed/24176112 输出的表达矩阵,如下所示:

请务必做出此题,很重要!

先给一下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)

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-02-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏增长技术

ffmpeg Documentation

| | | | | input | demuxer | encoded data | ...

571
来自专栏编程

Elasticsearch6.0 IKAnalysis分词使用

Elasticsearch 内置的分词器对中文不友好,会把中文分成单个字来进行全文检索,不能达到想要的结果,在全文检索及新词发展如此快的互联网时代,IK可以进行...

2396
来自专栏简书专栏

基于xgboost的风力发电机叶片结冰分类预测

xgboost中文叫做极致梯度提升模型,官方文档链接:https://xgboost.readthedocs.io/en/latest/tutorials/mo...

551
来自专栏简书专栏

基于tensorboard的模型训练过程可视化

2018年9月14日笔记 阅读本文的前提是已经阅读《基于tensorflow的一元二次方程回归预测》,文章链接:https://www.jianshu.com...

813
来自专栏zhangdd.com

ceph性能测试

该工具的语法为:rados bench -p <pool_name> <seconds> <write|seq|rand> -b <block size> -t...

1162
来自专栏Jack-Cui

Caffe学习笔记(四):使用pycaffe生成train.prototxt、test.prototxt文件

Python版本: Python2.7 运行平台: Ubuntu14.04 一、前言     了解到上一篇笔记的内容,就可以尝试自己编写python程序生...

5246
来自专栏州的先生

使用Django构建Python Restful Web服务:二、生成数据模型

1313
来自专栏生信技能树

得到一个物种所有基因的TSS(转录起始位点)区域的bed文件。

首先在UCSC的table browser 里面下载下面这个文件: ? 可以看到我这里选择的mm10的refseq系统的所有基因,共有29037个不同的tss,...

3008
来自专栏软件测试经验与教训

Python学习笔记(17)- os\\os.path 操作文件

3166
来自专栏PaddlePaddle

PaddlePaddle发布v0.10.0版

我们非常高兴发布了PaddlePaddle V0.10.0版,并开放了新的Python API。 之前在v0.9.0版,完成一个训练或预测任务至少需要两份pyt...

3507

扫码关注云+社区