从NCBI下载scRNA-seq数据后，如何用Python处理大体积fastq文件？

一、先做前置预处理（砍文件体积是关键） #

35GB的原始FASTQ直接用Python硬读会把内存撑爆，先通过命令行工具做质控过滤，大幅缩小文件大小：

• 用fastp做adapter修剪、低质量reads过滤，示例命令：

  fastp -i input.fastq -o cleaned.fastq -q 20 -u 30 --length_required 50
  

  这个命令会过滤掉质量值低于20的碱基占比超30%、长度不足50的reads，处理后文件体积能减少30%-60%。

二、Python处理大FASTQ的核心：流式/分块读取 #

绝对不要一次性把整个文件读进内存，用迭代器逐批处理：

1. 原生Python逐块处理 #

利用FASTQ每4行对应一条read的格式，批量读取处理：

def process_large_fastq(fastq_path, batch_size=10000):
    """批量处理FASTQ，每次处理batch_size条reads"""
    with open(fastq_path, 'r') as f:
        while True:
            batch = []
            # 读取一个批次的reads
            for _ in range(batch_size):
                lines = [f.readline().strip() for _ in range(4)]
                if not all(lines):
                    break
                batch.append({
                    'id': lines[0].lstrip('@'),
                    'seq': lines[1],
                    'qual': lines[3]
                })
            if not batch:
                break
            # 在这里写你的处理逻辑，比如统计GC含量、过滤reads
            process_batch(batch)

def process_batch(batch):
    # 示例：统计当前批次的平均GC含量
    gc_counts = sum((seq.count('G') + seq.count('C'))/len(seq) for seq in [b['seq'] for b in batch])
    avg_gc = gc_counts / len(batch)
    print(f"Batch average GC: {avg_gc:.2f}%")

2. 用pysam库简化流式读取 #

pysam专门处理测序数据，支持流式遍历FASTQ，内存占用极低：

import pysam

# 流式遍历每条read
with pysam.FastxFile('cleaned.fastq') as fh:
    for read in fh:
        # 获取read的核心信息
        read_id = read.name
        sequence = read.sequence
        quality = read.quality
        # 示例：过滤掉N碱基占比超10%的reads
        if sequence.count('N')/len(sequence) < 0.1:
            # 可以写入新的FASTQ文件，或者做其他处理
            pass

3. Dask分块做结构化分析 #

如果需要对FASTQ做统计分析（比如统计长度分布、碱基频率），用Dask分块加载成DataFrame，避免内存溢出：

import dask.dataframe as dd
import pandas as pd

def parse_fastq_chunk(chunk):
    """把分块的FASTQ行转换成DataFrame"""
    records = []
    for i in range(0, len(chunk), 4):
        if i+3 >= len(chunk):
            break
        records.append({
            'read_id': chunk.iloc[i].strip().lstrip('@'),
            'sequence': chunk.iloc[i+1].strip(),
            'quality': chunk.iloc[i+3].strip()
        })
    return pd.DataFrame(records)

# 按1GB分块读取FASTQ
dask_df = dd.read_csv(
    'cleaned.fastq',
    blocksize='1GB',
    header=None,
    sep='\n'
).map_partitions(parse_fastq_chunk)

# 计算全局平均read长度
avg_read_length = dask_df['sequence'].str.len().mean().compute()
print(f"Average read length: {avg_read_length:.1f}")

三、下游分析的高效路径 #

原始FASTQ处理完后，建议先做比对和定量：

• 用cellranger（针对10x数据）或STAR+featureCounts完成比对，得到稀疏基因表达矩阵（.mtx格式）
• 再用Python的scanpy或seurat-python分析，这些库默认用稀疏矩阵存储，能高效处理大样本量的数据，比如：

  import scanpy as sc
  
  # 读取稀疏矩阵
  adata = sc.read_10x_mtx('path_to_matrix')
  # 常规单细胞分析流程：过滤、归一化、降维
  sc.pp.filter_cells(adata, min_genes=200)
  sc.pp.normalize_total(adata)
  sc.pp.log1p(adata)
  sc.pp.pca(adata, n_comps=50)
  

四、硬件与性能优化 #

• 如果本地内存小于64GB，优先用云服务器的高内存实例（比如128GB内存），避免卡顿
• 用多进程加速批量处理，比如用concurrent.futures.ProcessPoolExecutor并行处理不同批次的reads

内容的提问来源于stack exchange，提问作者Derek