基因存在或缺失判断依据

3K文章判断依据，CDS>95%,body>85%

怎么判断测序样本的reads与参考基因组比对，每一个基因的CDS区域与基因主体区域，如果CDS区域大于95%且基因主体大于85%的判断为基因存在于样本中，否则该基因不存在？

1. reads与基因组进行比对。
2. 根据比对结果，提取reads在基因组上的位置信息。
3. 根据基因注释文件，确定CDS区域和基因主体的范围。
4. 计算reads在CDS区域的比例，判断是否大于95%。
5. 计算reads在基因主体的比例，判断是否大于85%。
6. 遍历比对结果，统计每个基因的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