3K文章判斷依據(jù)，CDS>95%,body>85%

怎么判斷測序樣本的reads與參考基因組比對，每一個基因的CDS區(qū)域與基因主體區(qū)域，如果CDS區(qū)域大于95%且基因主體大于85%的判斷為基因存在于樣本中，否則該基因不存在？

1. reads與基因組進行比對。
2. 根據(jù)比對結果，提取reads在基因組上的位置信息。
3. 根據(jù)基因注釋文件，確定CDS區(qū)域和基因主體的范圍。
4. 計算reads在CDS區(qū)域的比例，判斷是否大于95%。
5. 計算reads在基因主體的比例，判斷是否大于85%。
6. 遍歷比對結果，統(tǒng)計每個基因的CDS區(qū)域和基因主體區(qū)域的比對情況。

###下面是一個示例代碼，用于判斷每個基因的CDS區(qū)域和基因主體區(qū)域的比對情況，并判斷基因是否存在于樣本中：

```python
# 導入所需的庫
import pysam

# 輸入文件路徑
genome_file = "genome.fasta"  # 基因組文件
reads_file = "reads.fastq"  # reads文件
annotation_file = "annotation.gtf"  # 基因組注釋文件

# 定義變量
genes = {}  # 存儲每個基因的CDS區(qū)域和基因主體區(qū)域的比對情況

# 讀取基因組注釋文件，提取CDS區(qū)域和基因主體區(qū)域的位置信息
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文件，比對到基因組上并計數(shù)
with pysam.FastxFile(reads_file) as reads:
    for read in reads:
        seq = read.sequence
        # 進行比對，這里假設使用Bowtie進行比對，需提前安裝并設置好環(huán)境
        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區(qū)域和基因主體區(qū)域的比例，并判斷基因是否存在于樣本中
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區(qū)域比例：{:.2f}%".format(gene_info['cds_ratio']))
    print("基因主體比例：{:.2f}%".format(gene_info['gene_body_ratio']))
    print("是否存在于樣本中：{}".format(gene_info['exists_in_sample']))
    print()

來自知網(wǎng)

>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