功能:整理tRNAscan-SE输出的文件(-m),得到氨基酸和相应的tRNA的个数
为什么有这个需求:
高质量MAG还应编码23S、16S和5S rRNA基因,以及20种可能氨基酸中至少18种的tRNA
[文章:Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea]
Step1. 安装tRNAscan-SE
直接在Linux服务器用conda 安装
conda create -n trnascan
conda activate trnascan
conda install -c bioconda trnascan-se
题外话:我之前尝试了用下载源码安装——但这样以下面步骤这样安会有一个大问题,就是inferal没有同步装在所需目录,所以调用时就会报错。直接用conda安装就很方便地一切OK。
Step2. 运行 tRNAscan-SE
conda activate trnascan
# 这个test.fna是输入文件
tRNAscan-SE test.fna -B -o tRNA.out -f tRNA.ss -m tRNA.stats
得到的tRNA.stats,62行以后就是这样:
Step3. 数据整理(用R语言)
#######对于一个示例样品#########
# 读取文件内容
content <- readLines("tRNA.stats")
# 提取 Isotype / Anticodon Counts 以后的部分
start_line <- grep("Isotype / Anticodon Counts:", content) #得到这个定位的行数
extracted_content <- content[(start_line + 2):(length(content)-1)]
# 提取第一个制表符前的部分
extracted_content3 <- sub("\t.*", "", extracted_content)
# 提取第一个冒号前的部分(也即氨基酸的种类)
extracted_content2 <- sub(":.*", "", extracted_content3)
# 提取第一个制表符前,第一个冒号后的部分(也即比对到的个数)
extracted_content4 <- sub(".*:", "", extracted_content3)
df <- data.frame(Isotype = extracted_content2, count=extracted_content4,stringsAsFactors = FALSE)
输出的df表如下: