时间戳

首页 > 解决方案 > 从列表中获取计数数据后创建二进制矩阵文件

r - 从列表中获取计数数据后创建二进制矩阵文件

问题描述

在一个目录下，我有很多文本文件，文件的名称和格式如下：

SRS011061.txt
--------------------
contig SRS011061_idxstats.txt
BGC0000505 1
BGC0000505 1

SRS011090.txt
---------------------
contig SRS011090_idxstats.txt
BGC0000509 0
BGC0000509 1

SRS011271.txt
--------------------
contig SRS011271_idxstats.txt
BGC00001105 0
BGC00001105 0

从这些文件中，我想要两种类型的结果：

  1. 在每个文件中，当我们考虑不同的基因存在百分比（1/0 是存在和不存在）时，存在多少个基因（BGC**）。

对于这一步，我使用：

setwd("~/Desktop/test")
filenames <- list.files(full.names=F, pattern=".txt")
output <-lapply(filenames,function(i){
  t<-read.csv(i, header=T, check.names = F, sep = " ")
  t$gene_count<-1
  t[,2][t[,2]>0]<-1
  pre_abs<-aggregate(. ~ contig, t, sum)
  colnames(pre_abs)<-c("BGC_Accession","Gene_presence", "Gene_count")
  pre_abs$Percentage<-(pre_abs$Gene_presence/pre_abs$Gene_count)*100
  gene_pre_100_percent<-length(pre_abs$Percentage[pre_abs$Percentage>=100])
  gene_pre_20_percent<-length(pre_abs$Percentage[pre_abs$Percentage>=20])
  data.frame(Combinations=i,gene_pre_100_percent=gene_pre_100_percent,gene_pre_20_percent=gene_pre_20_percent)
})
Step2_TP_FP<-do.call(rbind,output)
Step2_TP_FP[,1] <- data.frame(gsub(".txt.*$", "", Step2_TP_FP[,1]))

在这里，我得到的结果是，在 100% 和 20% 的基因完成情况下，每个文件中存在多少个基因。

Combinations gene_pre_100_percent gene_pre_20_percent
SRS011061 1 1
SRS011090 0 1
SRS011271 0 0
  1. 对于相同的百分比标准，我想制作一个二进制存在缺席矩阵。对于 100% 标准，矩阵应如下所示：
SRS011061 SRS011090 SRS011271
BGC0000505 1 0 0
BGC0000509 0 0 0
BGC00001105 0 0 0

对于 20% 的标准，它应该是：

SRS011061 SRS011090 SRS011271
BGC0000505 1 0 0
BGC0000509 0 1 0
BGC00001105 0 0 0

得到step1结果后如何创建这两个矩阵文件？谢谢！

标签： r

解决方案

我认为，如果您保留所有信息并将存在率存储为百分比，而不是硬编码 20% 或 80%，然后一次性制作矩阵，这可能会更容易。例如对于一个文件，我们可以这样做：

output = lapply(fl,function(i){
   x = read.table(i,header=TRUE)
   cbind(aggregate(x[,2],list(gene=x[,1]),mean),
         file=sub(".txt","",i)
         )
})

output = do.call(rbind,output)

         gene   x      file
1  BGC0000505 1.0 SRS011061
2  BGC0000509 0.5 SRS011090
3 BGC00001105 0.0 SRS011271

现在只需设置获取矩阵的截止值，然后将其更改为 50%、60%，无需再次读取文件：

matrixfunc = function(da,perc){
  table(da$gene[da$x>perc],da$file[da$x>perc])
}

matrixfunc(output,0.8)
             
              SRS011061 SRS011090 SRS011271
  BGC0000505          1         0         0
  BGC0000509          0         0         0
  BGC00001105         0         0         0

matrixfunc(output,0.2)
             
              SRS011061 SRS011090 SRS011271
  BGC0000505          1         0         0
  BGC0000509          0         1         0
  BGC00001105         0         0         0

推荐阅读