prokka注释tsv文件统计
2021-06-29 本文已影响0人
灵木er
Prokka: rapid prokaryotic genome annotation
| .tsv | Tab-separated file of all features: locus_tag,ftype,len_bp,gene,EC_number,COG,produc |
|---|
我们有很多这样的prokka注释的tsv文件，怎么按照文件名给其分类汇总呢？单个prokka注释文件如下表所示：
| locus_tag | ftype | length_bp | gene | EC_number | COG | product |
|---|---|---|---|---|---|---|
| EFEEIJML_00001 | CDS | 447 | truB | 5.4.99.25 | tRNA pseudouridine synthase B | |
| EFEEIJML_00002 | CDS | 534 | hypothetical protein | |||
| EFEEIJML_00003 | CDS | 714 | hypothetical protein | |||
| EFEEIJML_00004 | CDS | 315 | hypothetical protein | |||
| EFEEIJML_00005 | CDS | 639 | hypothetical protein | |||
| EFEEIJML_00006 | CDS | 627 | rnhB | 3.1.26.4 | COG0164 | Ribonuclease HII |
| EFEEIJML_00007 | CDS | 708 | hypothetical protein | |||
| EFEEIJML_00008 | CDS | 963 | ndhF | 1.17.1.5 | Nicotinate dehydrogenase FAD-subunit | |
| EFEEIJML_00009 | CDS | 1551 | 2.1.1.156 | Glycine/sarcosine/dimethylglycine N-methyltransferase |
最终输出表如下：
| file_ID | Total_num | EC_num | COG_num | GENE_num | CDS_num | rRNA_num | tRNA_num | tmRNA_num |
|---|---|---|---|---|---|---|---|---|
| file.1 | 3131 | 1006 | 1187 | 1478 | 3089 | 3 | 38 | 1 |
| file.10 | 1639 | 651 | 744 | 952 | 1617 | 0 | 22 | 0 |
| file.100 | 3600 | 1354 | 1626 | 1973 | 3550 | 2 | 48 | 0 |
| file.101 | 2723 | 1050 | 1421 | 1776 | 2684 | 6 | 32 | 1 |
| file.102 | 3524 | 1146 | 1224 | 1586 | 3491 | 0 | 32 | 1 |
| file.103 | 2063 | 739 | 991 | 1136 | 2020 | 0 | 42 | 1 |
| file.104 | 2478 | 874 | 1049 | 1353 | 2438 | 1 | 39 | 0 |
| file.105 | 2955 | 1049 | 1352 | 1710 | 2926 | 1 | 27 | 1 |
| file.106 | 2215 | 732 | 979 | 1197 | 2184 | 0 | 30 | 1 |
excel的分类汇总应该可以吧，不太熟。这里使用R语言来实现，涉及到批量读取文件并合并，分类汇总，主要就是这些。
代码如下：
pacman::p_load(purrr, stringr, data.table)
files = list.files("prokka_out", pattern="tsv", full.names=T, recursive=T)
df1 = map_dfr(set_names(files), fread, .id="file")
cog = df1[,length(COG)- sum(COG==""), by=(file)] %>% .[order(file)]
ec = df1[,length(EC_number)- sum(EC_number==""), by=(file)] %>% .[order(file)]
gene = df1[,length(gene)- sum(gene==""), by=(file)] %>% .[order(file)]
total = df1[, .N, by=.(file)] %>% .[order(file)]
ftype = df1[, .N, by=.(file, ftype)] %>% dcast(file ~ ftype, value="N", fill=0) %>% .[order(file)]
fileid = df1[, str_extract(file, "file\\.\\d*"), by=.(file)] %>% .[order(file)]
data1 = cfiled(fileid,total,ec,gene,cog,ftype)
data1 = data1[, !"file"]
names(data1) = c("file_ID", "Total_num", "EC_num", "GENE_num", "COG_num", "CDS_num", "rRNA_num", "tRNA_num", "tmRNA_num")
fwrite(data1, "prokka_table.tsv", sep="\t")
代码只有十几行，逻辑也很清楚，没有用到循环，
第一行：加载需要的包，purrr提供map_系列函数及管道%>%，stringr处理字符串，这里是提取字符串，data.table提供高性能读写及计算；
第二行：构建待读取文件名；
第三行：使用fread读文件，map_dfr迭代并合并；
第四-九行：按照file分类汇总，order排序是为了方便直接合并，不排序可以使用连接，dcast实现长变宽；
第十行：直接合并；
第十一行：剔除不需要的行，连接的时候分类汇总的group重复了多次，反选剔除，也可以在合并的时候去掉第一列；
第十二行：表头重命名；
第十三行：写入文件。
这里主要是使用data.table包的操作处理的，也可以使用read_tsv函数读写，使用tidyverse处理。文件不大使用那个都行，时间差别不大，文件比较大则建议使用data.table包，快！！！
参考：张敬信. R 语言编程 — 基于 tidyverse[M]. .
