3K文章判断依据，CDS>95%,body>85%
怎么判断测序样本的reads与参考基因组比对，每一个基因的CDS区域与基因主体区域，如果CDS区域大于95%且基因主体大于85%的判断为基因存在于样本中，否则该基因不存在？
- reads与基因组进行比对。
- 根据比对结果，提取reads在基因组上的位置信息。
- 根据基因注释文件，确定CDS区域和基因主体的范围。
- 计算reads在CDS区域的比例，判断是否大于95%。
- 计算reads在基因主体的比例，判断是否大于85%。
- 遍历比对结果，统计每个基因的CDS区域和基因主体区域的比对情况。
###下面是一个示例代码，用于判断每个基因的CDS区域和基因主体区域的比对情况，并判断基因是否存在于样本中：
```python
# 导入所需的库
import pysam
# 输入文件路径
genome_file = "genome.fasta" # 基因组文件
reads_file = "reads.fastq" # reads文件
annotation_file = "annotation.gtf" # 基因组注释文件
# 定义变量
genes = {} # 存储每个基因的CDS区域和基因主体区域的比对情况
# 读取基因组注释文件，提取CDS区域和基因主体区域的位置信息
with open(annotation_file, 'r') as f:
for line in f:
if not line.startswith('#'):
fields = line.strip().split('\t')
feature_type = fields[2]
if feature_type == 'gene':
gene_id = fields[8].split(';')[0].split(' ')[1].replace('"', '')
genes[gene_id] = {'cds_count': 0, 'gene_body_count': 0, 'cds_ratio': 0, 'gene_body_ratio': 0}
elif feature_type == 'CDS':
start = int(fields[3])
end = int(fields[4])
gene_id = fields[8].split(';')[0].split(' ')[1].replace('"', '')
if gene_id in genes:
genes[gene_id]['cds_regions'] = genes[gene_id].get('cds_regions', []) + [(start, end)]
elif feature_type == 'exon':
start = int(fields[3])
end = int(fields[4])
gene_id = fields[8].split(';')[0].split(' ')[1].replace('"', '')
if gene_id in genes:
genes[gene_id]['gene_body_regions'] = genes[gene_id].get('gene_body_regions', []) + [(start, end)]
# 遍历reads文件，比对到基因组上并计数
with pysam.FastxFile(reads_file) as reads:
for read in reads:
seq = read.sequence
# 进行比对，这里假设使用Bowtie进行比对，需提前安装并设置好环境
sam_output = pysam.view("-S", "-q", "10", genome_file, "-", "-o", "temp.sam", "-")
with open("temp.sam", 'r') as samfile:
for line in samfile:
fields = line.strip().split('\t')
gene_id = fields[2] # 假设基因信息在SAM文件的第三列
if gene_id in genes:
genes[gene_id]['cds_count'] += 1 if any(start <= int(fields[3]) <= end for start, end in genes[gene_id]['cds_regions']) else 0
genes[gene_id]['gene_body_count'] += 1 if any(start <= int(fields[3]) <= end for start, end in genes[gene_id]['gene_body_regions']) else 0
# 计算每个基因的CDS区域和基因主体区域的比例，并判断基因是否存在于样本中
for gene_id, gene_info in genes.items():
cds_ratio = gene_info['cds_count'] / total_count * 100
gene_body_ratio = gene_info['gene_body_count'] / total_count * 100
gene_info['cds_ratio'] = cds_ratio
gene_info['gene_body_ratio'] = gene_body_ratio
if cds_ratio > 95 and gene_body_ratio > 85:
gene_info['exists_in_sample'] = True
else:
gene_info['exists_in_sample'] = False
# 打印结果
for gene_id, gene_info in genes.items():
print("基因ID：{}".format(gene_id))
print("CDS区域比例：{:.2f}%".format(gene_info['cds_ratio']))
print("基因主体比例：{:.2f}%".format(gene_info['gene_body_ratio']))
print("是否存在于样本中：{}".format(gene_info['exists_in_sample']))
print()
来自知网
>80%外显子
使用SGSGeneLoss软件
参考：
Wang, W., Mauleon, R., Hu, Z. et al. Genomic variation in 3,010 diverse accessions of Asian cultivated rice. Nature 557, 43–49 (2018). https://doi.org/10.1038/s41586-018-0063-9
Gao, L., Gonda, I., Sun, H. et al. The tomato pan-genome uncovers new genes and a rare allele regulating fruit flavor. Nat Genet 51, 1044–1051 (2019). https://doi.org/10.1038/s41588-019-0410-2
