二代测序和三代测序的区别 然后基于imargi的测序结果要进行处理时,些bash代码的时候应该有什么...
생성일: 2026년 3월 4일
생성일: 2026년 3월 4일
二代测序和三代测序的区别 然后基于imargi的测序结果要进行处理时,些bash代码的时候应该有什么样的区别
下面分两块说:二代(NGS/短读长) vs 三代(TGS/长读长)的核心区别;以及iMARGI(你写成 imargi,我按 iMARGI 理解)数据在写 bash 处理时,和普通二代测序有哪些“必做的特殊处理点”。
读长(read length)
单分子 & PCR 偏好
下游软件生态(直接影响 bash 写法)
fastqc/multiqc -> cutadapt/fastp -> bwa mem/bowtie2/STAR -> samtools/picard -> (varcall/quant/peak...)basecalling(ONT guppy/dorado 或 PacBio) -> NanoPlot/pycoQC -> minimap2/pbmm2 -> samtools -> (polish/iso-seq/SV...)一句话:二代 bash 更多围绕“短读长的清洗+比对+去重复/定量”;三代 bash 更多围绕“basecall/长读长过滤+长读长比对+抛光/结构变异/全长转录本”。
iMARGI 的核心是把一次 RNA–DNA 相互作用变成一个RNA–linker–DNA 的嵌合分子,再用paired-end去读出(通常 Illumina PE)。(eScholarship)
因此它虽然是二代 PE 数据,但**清洗、比对参数、以及“把 BAM 解析成 RNA–DNA pairs”**都和常规 WGS/RNA-seq 不一样。
4DN/iMARGI 管线里明确写了两点:(data.4dnucleome.org)
seqtk trimfq -b 2这两点会直接改变你 bash 里的“trim/filter”逻辑:
4DN 管线给的关键命令:(data.4dnucleome.org)
bashbwa mem -t <nthreads> -SP5M <genome_index> <fastq1> <fastq2>
这些参数并不是“随便写写”:
-SP:让结果更接近“分别比对两端”但仍保留 PE 格式(避免把两端强行当成同一基因组片段去牵引纠错)-5:对嵌合比对(chimeric alignments)时,优先把 5' 端作为主比对(后续解析更一致)(data.4dnucleome.org)普通二代(比如 WGS)你可能就 bwa mem 默认参数 + markdup;但 iMARGI 建议严格按这个来。
4DN 总结 iMARGI 管线的核心组件就是:clean & align -> parse 成 pairs -> merge/aggregate(甚至生成 cool)。(data.4dnucleome.org)
这一步是 iMARGI 和普通二代最大的本质差别:
另外一个容易写错的点:Zhong lab 提醒 RNA 端(R1)是反向链特异(reverse strand specific),你后面做基因注释/方向时要把链翻转;DNA 端不是链特异。(GitHub)
下面是骨架示例(你可以替换样本名/路径/线程数)。其中 CT 过滤用 Zhong 方法(更严格),R1 trim用 4DN/seqtk 思路。
bash#!/usr/bin/env bash set -euo pipefail THREADS=16 REF=./ref/GRCh38.fa BWAIDX=./ref/bwa_index/GRCh38 R1=./raw/sample_R1.fastq.gz R2=./raw/sample_R2.fastq.gz mkdir -p clean mapping pairs # 1) iMARGI-clean: R1 去掉 5' 端 2bp 随机碱基(4DN 推荐) zcat "$R1" | seqtk trimfq -b 2 - | gzip > clean/sample.clean.R1.fastq.gz # 2) iMARGI-clean: 过滤 R2 前两碱基不是 "CT" 的 read pair(Zhong 方法) # 注:真实生产建议用 Zhong 提供的 ct_clean_pair.sh(含并行/配对修复)([GitHub](https://github.com/Zhong-Lab-UCSD/iMARGI_methods)) paste \ <(zcat clean/sample.clean.R1.fastq.gz | paste - - - -) \ <(zcat "$R2" | paste - - - -) | \ awk '{ # R2 的序列行在第 6 列(因为每端 4 行被 paste 成 4 列) r2seq=$6 if (substr(r2seq,1,2)=="CT") print }' | \ tee >(cut -f1-4 | tr '\t' '\n' | gzip > clean/sample.clean.CT.R1.fastq.gz) \ >(cut -f5-8 | tr '\t' '\n' | gzip > clean/sample.clean.CT.R2.fastq.gz) \ > /dev/null # 3) mapping: bwa mem -SP5M(4DN 推荐)([data.4dnucleome.org](https://data.4dnucleome.org/resources/data-analysis/imargi-pipeline)) bwa index -p "$BWAIDX" "$REF" 2>/dev/null || true bwa mem -t "$THREADS" -SP5M "$BWAIDX" \ clean/sample.clean.CT.R1.fastq.gz \ clean/sample.clean.CT.R2.fastq.gz | \ samtools view -Shb - | \ samtools sort -@ "$THREADS" -o mapping/sample.bam samtools index mapping/sample.bam # 4) parse -> pairs: # 这里通常调用 iMARGI 官方/4DN 的解析脚本或 docker pipeline(会输出 pairs / cool 等)([data.4dnucleome.org](https://data.4dnucleome.org/resources/data-analysis/imargi-pipeline)) # 例如(示意): # imargi_parse ... -i mapping/sample.bam -o pairs/sample.pairs.gz
-SP5M (data.4dnucleome.org)如果你愿意把你现在的目录结构、read 长度(2×100 还是 2×150)、物种(hg38/mm10)、以及你准备输出的是 pairs 还是 cool发我(不用发数据本身,只要这些参数),我可以把上面的骨架改成一份你机器上能直接跑的“完整脚本模板”(含日志、并行、失败退出、multiqc 汇总)。