刘小泽写于2018.8.20
利用htseq-count、featureCounts等定量软件生成的counts计数文件都是单独的，后续进行统计需要将它们汇集到一起
实例测试（单行perl结合R语言）
现在有三个count文件，SRR3589956.count、SRR3589957.count、SRR3589958.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命令
得到的结果
第三步 探索
- 统计各个样本count值（最值、中位数）
做个统计
-
按照基因来统计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!
网友评论