这里使用111份文章里构建的pan-genome，实际注释到的基因条目有154010个
genomecov 示意图
文章中使用
genomecov -bga -split参数统计覆盖度，实际上是统计reads的叠加情况，试了一下，并不能直接统计基因区间的具体覆盖情况
二代需要>85%，前面可能理解错了
这里使用bedtools coverage,优势是可以直接计算出目标区域的平均深度。
首先从注释文件中得到gene/CDS区间的bed位置文件
gene区间，分为3列，tab分割
一个gene的CDS区间实际上过于抽象，
统计窗口内的平均覆盖深度
bedtools coverage -a gene.bed -b bams/DL651.bam >cov/DL651.cov
实际>95%覆盖度的条目有53074；>85%有75177个
每一列分别代表染色体、起止位点、reads数、碱基数、区间大小/基因长度、覆盖率
这样就可以大致确定每个基因的PAV情况了。
因为每一行的顺序是不变的，所有提取第七列方便后面合并，当然也可以一步完成
cat fqid|while read id
do
cut -f 7 ${id}.cov >tsv/$id.tsv
done
从只含有gene的注释提取id信息
###tq1.py
import re
# 读取文件
with open('gene', 'r') as file:
lines = file.readlines()
# 提取第九列ID后的字符
results = []
for line in lines:
columns = line.strip().split('\t') # 假设文件是以制表符分隔的
if len(columns) >= 9:
match = re.search(r'ID=(.*?);', columns[8]) # 第九列是索引为8
if match:
results.append(match.group(1))
# 输出结果
for result in results:
print(result)
###python tq1.py >geneid
整合完之后大概是这样子
替换成0，1
# 读取文件
data <- read.table("allcov.tsv", header = TRUE, sep="\t")
# 从第二列第二行开始，将符合条件的值替换
data[1:nrow(data), 2:ncol(data)] <- ifelse(data[1:nrow(data), 2:ncol(data)] > 0.85, 1, 0)
# 保存修改后的数据框
write.table(data, "allcov_modified.tsv", sep="\t", quote=FALSE, row.names=FALSE)
替换后的结果
