表达矩阵的合并

刘小泽写于2018.8.20
利用htseq-count、featureCounts等定量软件生成的counts计数文件都是单独的,后续进行统计需要将它们汇集到一起

实例测试(单行perl结合R语言)

现在有三个count文件,SRR3589956.countSRR3589957.countSRR3589958.count

《表达矩阵的合并》 每一个count文件长这样

第一步 初步合并各个计数文件

先将两个文件合并,第一列是样本名,第二列是基因名,第三列是count结果

perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count >matrix.count

《表达矩阵的合并》 三个先合并到一个

第二步 R重塑合并的计数矩阵

a=read.csv('matrix.count',header=F,sep="\t")
colnames(a)=c('sample','gene','count')
library(reshape2)
counts=dcast(a,formula=gene~sample)
write.table(counts, file="join.count",sep="\t",quote=FALSE,row.names=FALSE)

可以直接在命令行运行(前提是安装了R)

《表达矩阵的合并》 R命令
《表达矩阵的合并》 得到的结果

第三步 探索

  1. 统计各个样本count值(最值、中位数)

《表达矩阵的合并》 做个统计

  1. 按照基因来统计counts的平均数

    install.packages("dplyr")
    library(dplyr)
    gene_mean = group_by(a,gene)%>%summarize(mean=mean(count))
    #另外也可以按照sample来统计
    #当然,除了mean平均数,还可以求中值median,最值max、min
    

《表达矩阵的合并》 按基因来统计count

(没有实际数据)自己生成测试数据

目的:生成a-g.txt 7个文件,每个文件中第一列为基因名,从gene_1到gene_99,第二列是表达量,1000以内随机整数

#新建测试文件test.sh
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
# -e后面紧跟着引号里面的字符串是要执行的命令;
# -l输出换行
#想要输出gene_1这样类型的,gene_后面的数字用一个变量$_代替,这个变量相当于占一个位置,它的赋值是在后面foreach 1..99,表示 $_从1到99
#加一个\t表示一个tab键,空4格
#rand生成随机数,int限制整数
#注意到"gene_$_\t"与后面生成整数函数之间有一个.【不加.就会报错,它的作用应该是分离字符串和函数】
#perl脚本调用test.sh
perl -e 'system("bash test.sh > $_.txt") foreach a..g'
#结果就生成了想要的测试文件

开始统计:

#先合并
perl -lne 'if ($ARGV=~/(.*).txt/){print $1\t$_}' *.txt > matrix.count
#进入R studio
#然后重复上面👆第二步

欢迎关注我们的公众号~_~  
我们是两个农转生信的小硕,打造生信星球,想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

《表达矩阵的合并》 Welcome to our bioinfoplanet!

    原文作者:刘小泽
    原文地址: https://www.jianshu.com/p/853f3706eaa8
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞