不管是重测序或者是有参转录组，很重要的一环就是将测序数据比对到参考基因组上。一般来说，常用比对软件都可以正常运行结束，但是当碰到超大基因组的时候难免会碰到一些问题，这些问题有的时候是参考基因组索引建立方面的，有的是比对途中的，有的是比对结果异常的。

关于参考基因组索引建立：通常来说都是有于内存不足引起的，考虑租用更大的服务器、更换到更大内存的节点可以得到解决。

比对过程中少有问题，特别是重测序，由于需求量大，所以软件优化更新都挺及时的，一般可以通过软件升级/更换更优软件解决问题。

而比对结果异常通常大家可能碰到以下方面的问题：

1. 单条染色体长度超出，无法构建bai索引

我们常规构建bam文件的索引（bai）使用：

samtools index a.bam

此时需要增加 -c参数，构建 csi的索引：

samtools index -c a.bam

2. 比对文件中存在负坐标

这个问题并不多见，通常出现在超大基因组上，单条序列长度超过2.14G才会导致负坐标出现。

因为bam/sam存储位置信息时使用int32_t数据类型，因此它所能接纳的最大染色体长度为2^31 - 1 = 2,147,483,647 bp（约2.14G）。

此处我并没有找到更好的软件可以解决这个问题，因此只能对超长序列进行拆分。

比如在这里我需要把chr1拆掉。请注意，这里因为我们的序列长度超过2.14G所以无法使用bedtools的getfasta功能。大家可以自己写python脚本进行拆除，但是不要写成读入整个基因组之后才能进行拆除的模式。

samtools faidx ref.fasta 

根据上面的ref.fasta.fai决定怎么拆。更严谨的是统计fasta里面N的所在位置，然后在N的地方拆。次一等的是，大概选定一个所在位置之后，去gff3文件里面找一下，断点不要处于基因结构内部就可以。

seqkit subseq --chr chr1 -r 1:1200000000 ref.fasta > chr1_part.fa
seqkit subseq --chr chr1 -r 1200000001:2800000000 ref.fasta >> chr1_part.fa
seqkit grep -v -p "chr1" ref.fasta >> chr1_part.fa
samtools faidx chr1_part.fa

更改一下序列的名字，需要根据chr1_part.fa.fai去做一个bed，第一列是现在chr1_part.fa的染色体名字，第二列为0，第三列是终止位置，也就是chr1_part.fa.fai的第二列，第四列是你想让这条序列叫什么。比如我的bed第一行是chr1:1-1200000000 0 1200000000  chr1_part1：

bedtools getfasta -fi chr1_part.fa -bed bed  -nameOnly -fo ref.fa

根据比对软件的需求，对前面拆除的fa重新建立索引，然后完成比对就可以了。

如果是转录组的话，需要对gtf也进行一下切割，使用如下：

python gtf_cutter.py -g raw.gtf -c split.txt -o ref.gtf
#split.txt 就这一行：chr1    1200000000，如果有多个要切割的染色体就写多行，第二列是切割的位置。

gtf_cutter.py如下：

#!/usr/bin/env python3
"""
GTF染色体切割工具
专门用于根据切割位置处理GTF/GFF文件，调整染色体名称和特征坐标
"""
import argparse
import sys
import os
import gzip
from collections import defaultdict
def read_cut_positions(cut_file):
    """
    读取切割位置文件
    格式：染色体<TAB>切割位置
    例如：chr1<TAB>120000000
    """
    cut_dict = {}
    try:
        with open(cut_file, 'r') as f:
            for line_num, line in enumerate(f, 1):
                line = line.strip()
                if not line or line.startswith('#'):
                    continue
                parts = line.split('\t')
                if len(parts) < 2:
                    print(f"警告: 第{line_num}行格式不正确，跳过")
                    continue
                
                chrom = parts[0].strip()
                try:
                    cut_pos = int(parts[1].strip())
                    if cut_pos <= 0:
                        print(f"警告: 第{line_num}行切割位置必须为正整数，跳过")
                        continue
                    cut_dict[chrom] = cut_pos
                except ValueError:
                    print(f"警告: 第{line_num}行切割位置不是有效整数，跳过")
                    
        print(f"成功读取 {len(cut_dict)} 个切割位置")
        return cut_dict
        
    except FileNotFoundError:
        print(f"错误: 切割位置文件 '{cut_file}' 不存在")
        sys.exit(1)
    except Exception as e:
        print(f"读取切割位置文件时出错: {e}")
        sys.exit(1)
def process_gtf_file(gtf_file, cut_dict, output_gtf):
    """
    处理GTF/GFF文件，根据染色体切割调整坐标
    """
    try:
        # 自动检测是否为压缩文件
        if gtf_file.endswith('.gz'):
            open_func = gzip.open
            mode = 'rt'
        else:
            open_func = open
            mode = 'r'
        
        stats = defaultdict(int)  # 统计信息
        
        with open_func(gtf_file, mode) as fin, open(output_gtf, 'w') as fout:
            for line_num, line in enumerate(fin, 1):
                original_line = line
                line = line.strip()
                
                # 保留注释行原样输出
                if not line or line.startswith('#'):
                    fout.write(original_line)
                    stats['comments'] += 1
                    continue
                
                parts = line.split('\t')
                if len(parts) < 5:
                    fout.write(original_line)
                    stats['invalid'] += 1
                    continue
                
                chrom = parts[0]
                try:
                    start_pos = int(parts[3])
                    end_pos = int(parts[4])
                except ValueError:
                    # 如果坐标不是数字，保持原样
                    fout.write(original_line)
                    stats['invalid_coord'] += 1
                    continue
                
                stats['total_features'] += 1
                
                # 检查是否需要切割该染色体
                if chrom in cut_dict:
                    cut_pos = cut_dict[chrom]
                    stats['chromosomes_cut'] += 1
                    
                    # 特征完全在切割点之前
                    if end_pos <= cut_pos:
                        parts[0] = f"{chrom}_part1"
                        stats['features_before_cut'] += 1
                    
                    # 特征完全在切割点之后
                    elif start_pos > cut_pos:
                        parts[0] = f"{chrom}_part2"
                        parts[3] = str(start_pos - cut_pos)
                        parts[4] = str(end_pos - cut_pos)
                        stats['features_after_cut'] += 1
                    
                    # 特征跨越切割点，提供处理选项
                    elif start_pos <= cut_pos and end_pos > cut_pos:
                        # 选项1: 跳过跨越切割点的特征
                        # print(f"警告: 第{line_num}行特征跨越切割点，已跳过 (染色体 {chrom}, 位置 {start_pos}-{end_pos})")
                        # stats['features_spanning_cut'] += 1
                        # continue
                        
                        # 选项2: 分割特征（复杂，这里简单处理为分配到第一部分）
                        parts[0] = f"{chrom}_part1"
                        parts[4] = str(cut_pos)  # 将结束位置调整为切割点
                        stats['features_spanning_cut'] += 1
                        print(f"警告: 第{line_num}行特征跨越切割点，已截断 (染色体 {chrom}, 位置 {start_pos}-{end_pos})")
                    
                    # 重新组装行
                    new_line = '\t'.join(parts) + '\n'
                    fout.write(new_line)
                    stats['features_modified'] += 1
                    
                else:
                    # 不需要切割的染色体，原样输出
                    fout.write(original_line)
                    stats['chromosomes_unchanged'] += 1
        
        # 输出统计信息
        print(f"\n=== GTF处理统计 ===")
        print(f"总处理行数: {line_num}")
        print(f"注释行: {stats['comments']}")
        print(f"总特征数: {stats['total_features']}")
        print(f"需要切割的染色体数: {len(cut_dict)}")
        print(f"切割前特征数: {stats['features_before_cut']}")
        print(f"切割后特征数: {stats['features_after_cut']}")
        print(f"跨越切割点特征数: {stats['features_spanning_cut']}")
        print(f"修改的特征总数: {stats['features_modified']}")
        print(f"未修改的特征数: {stats['chromosomes_unchanged']}")
        print(f"无效行数: {stats['invalid'] + stats['invalid_coord']}")
        print(f"修改后的GTF文件已保存至: {output_gtf}")
        
    except Exception as e:
        print(f"处理GTF文件时出错: {e}")
        sys.exit(1)
def validate_gtf_file(gtf_file):
    """
    简单验证GTF文件格式
    """
    try:
        if gtf_file.endswith('.gz'):
            open_func = gzip.open
            mode = 'rt'
        else:
            open_func = open
            mode = 'r'
        
        with open_func(gtf_file, mode) as f:
            for i, line in enumerate(f):
                if i >= 10:  # 只检查前10行
                    break
                if line.strip() and not line.startswith('#'):
                    parts = line.split('\t')
                    if len(parts) >= 5:
                        print(f"GTF格式验证通过（找到 {len(parts)} 列）")
                        return True
        return True
    except Exception as e:
        print(f"GTF文件验证警告: {e}")
        return True
def main():
    parser = argparse.ArgumentParser(description='GTF染色体切割工具 - 专门处理GTF/GFF文件')
    parser.add_argument('-g', '--gtf', required=True, help='输入GTF/GFF文件路径（支持.gz格式）')
    parser.add_argument('-c', '--cuts', required=True, help='切割位置文件（TSV格式：染色体<TAB>位置）')
    parser.add_argument('-o', '--output', required=True, help='输出GTF文件路径')
    parser.add_argument('--validate', action='store_true', help='验证输入GTF文件格式')
    
    args = parser.parse_args()
    
    # 检查输入文件是否存在
    for file_path in [args.gtf, args.cuts]:
        if not os.path.exists(file_path):
            print(f"错误: 文件 '{file_path}' 不存在")
            sys.exit(1)
    
    print("=== GTF染色体切割工具开始运行 ===")
    
    # 可选：验证GTF文件格式
    if args.validate:
        validate_gtf_file(args.gtf)
    
    # 1. 读取切割位置
    cut_dict = read_cut_positions(args.cuts)
    
    if not cut_dict:
        print("警告: 未找到有效的切割位置，输出文件将与输入文件相同")
    
    # 2. 处理GTF文件
    process_gtf_file(args.gtf, cut_dict, args.output)
    
    print("=== GTF处理完成 ===")
if __name__ == "__main__":
    main()