bedtools coverage统计目标区间覆盖度

这里使用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)

替换后的结果