高效统计polyA中每个序列GCT碱基数量

发布于:2022-08-09 ⋅ 阅读:(371) ⋅ 点赞:(0)

 文件数据格式如下:

输出格式如下:

#读取csv文件
library(data.table)
library(stringr)
library(dplyr)
print(getwd())
start<-Sys.time()
polyAData <- fread(file="./data/D10.PolyA_FL.csv",encoding='UTF-8')
#输出行列数量值
cat(nrow(polyAData),ncol(polyAData))
#获取seq列数据
print("----------创建向量,统计polyA中GCT数据,后面合成数据框添加到原始数据框中----------")
Total = c()
G = c()
C = c()
T = c()
GCT = c()
print("-------------------apply,select")
#统计碱基相关数据
statiGCT=function(x){
  polyA_Total=nchar(x)
  g_count=str_count(x,pattern = "G")
  c_count=str_count(x,pattern = "C")
  t_count=str_count(x,pattern = "T")
  gct_count = g_count+c_count+t_count;
  Total<<-c(Total,polyA_Total)
  G<<-c(G,g_count)
  C<<-c(C,c_count)
  T<<-c(T,t_count)
  GCT<<-c(GCT,gct_count)

#取第二列数据,按列取,将该列中每行数据作为变量x传入statiGCT方法中
#apply(select(polyAData,2),2,function(x) statiGCT(x) )
apply(select(polyAData,2),2,function(x) statiGCT(x) )
#将向量组合成数据框
print("将向量组合成数据框")
GCTStati = cbind(Total,G,C,T,GCT)
write.csv(GCTStati,file = "./data/output/mydataTest5.csv",row.names = F)
end<-Sys.time()

runningtime<-end-start
print(cat("程序运行时间:",runningtime) )