戰(zhàn):7 個高頻疑問帶你跑出第一份 SNP 檢出)
Snippy 變異檢測入門實(shí)戰(zhàn)7 個高頻疑問帶你跑出第一份 SNP 檢出【免費(fèi)下載鏈接】snippy:scissors: :zap: Rapid haploid variant calling and core genome alignment項目地址: https://gitcode.com/gh_mirrors/sn/snippySnippy 變異檢測工具能把測序 reads 對齊到參考基因組快速找出 SNP 與插入缺失還能把多個樣本對齊到同一套參考上做核心基因組比較。這篇不按說明書順序講而是把新手最常問的 7 個問題挨個拆開每個問題都配先做什么、看到什么、下一步的可對照步驟。 問題一裝 Snippy 有哪三條路我該怎么選現(xiàn)實(shí)里沒有第四條官方默認(rèn)路線主流做法就三種Conda、Homebrew、源碼。它們之間的核心差異其實(shí)是**依賴由誰來管**的差異。先看這張自己設(shè)計的三選一對比表安裝路線一句話定位適合誰要留意的坑Conda一條命令把 Snippy 和全部依賴一起裝齊想省事、不想跟依賴?yán)p斗的人包版本可能比上游略舊Homebrew用 brew 統(tǒng)一管理和系統(tǒng)其他軟件一視同仁macOS 用戶Linux 可配 LinuxBrew依賴按 brew 規(guī)則解析個別工具可能要手動補(bǔ)源碼從倉庫拉最新代碼永遠(yuǎn)不過時追新功能、愿意折騰的人PATH 自己配依賴全部自己裝路線 Aconda 裝 Snippyconda install -c conda-forge -c bioconda -c defaults snippy敲下去之后終端會列出一串待安裝清單里面除了 snippy 本體還會出現(xiàn) bwa、freebayes、samtools 這些名字。進(jìn)度條走完、回到命令提示符就算成了。路線 Bbrew 一條命令brew install brewsci/bio/snippy適合平時就習(xí)慣用 brew 管理軟件的用戶。路線 C源碼永遠(yuǎn)最新git clone https://gitcode.com/gh_mirrors/sn/snippy.git export PATH$PWD/snippy/bin:$PATH第二行的路徑要替換成你實(shí)際克隆下來的位置。跑完snippy --help能彈出完整參數(shù)說明而不是command not found就是成功。帶走什么不想操心依賴就選 condamacOS 用戶優(yōu)先 brew非要最新功能才走源碼。 問題二裝好 Snippy 本體為什么還不能直接開工這是新手最容易踩的隱形坑。Snippy 的定位更像一個包工頭它負(fù)責(zé)排工序、派活真正砌墻的是它手下那批外部工具——bwa 負(fù)責(zé)比對、freebayes 負(fù)責(zé)變異識別、samtools/bcftools 處理 BAM 與 VCF、snpEff 負(fù)責(zé)注釋背后還跟著 bedtools、seqtk、samclip 等一票配角。缺了任何一個流程都會在某個步驟卡住或報錯。所以選路線時你真正該關(guān)注的問題不是Snippy 本體裝沒裝上而是**它的依賴能不能被一并解決**。這也是 conda 路線對新手最友好的原因——bioconda 會把整條依賴鏈一次補(bǔ)齊。帶走什么把裝 Snippy理解為裝一個帶全套幫手的工具箱心態(tài)就對了。? 問題三裝完怎么用兩個命令驗收環(huán)境別急著上真實(shí)數(shù)據(jù)先花兩分鐘驗貨。第一步確認(rèn)版本snippy --version當(dāng)前倉庫對應(yīng)的版本號是snippy 5.0.0-dev。能打印出這串字符說明本體就位報錯的話回頭查 PATH。第二步檢查全部依賴snippy --check它會逐個探測 bwa、minimap2、samtools、bcftools、bedtools、freebayes、snpEff 等組件每個可用項旁邊顯示 OK 或?qū)?yīng)的版本信息。哪個標(biāo)成缺失就針對哪個補(bǔ)conda install -c bioconda samtools bcftools bwa freebayes snpeff samclip seqtk補(bǔ)完再跑一遍snippy --check直到全部通過再往下走。帶走什么這兩條命令就是你的開工許可比任何安裝日志都可信。 問題四我的數(shù)據(jù)到底適不適合喂給 Snippy先劃適用范圍Snippy 面向單倍體基因組。細(xì)菌、病毒、質(zhì)粒、線粒體這些每個位點(diǎn)只有一份拷貝的樣本都是它的主場。反過來人類的二倍體數(shù)據(jù)不適合它——二倍體要分辨雜合位點(diǎn)那不是 Snippy 的設(shè)計目標(biāo)。再看輸入形態(tài)Snippy 其實(shí)相當(dāng)好說話雙端 FASTQ用--R1/--R2傳進(jìn)去單端 FASTQ只給--R1也能跑拼裝好的 contigs走--ctgs后面的場景部分會專門講參考基因組可以是 FASTA 或 GENBANKreads 支持 gz 壓縮。把這幾條記牢實(shí)戰(zhàn)時就不會在喂數(shù)據(jù)這一步卡殼。帶走什么單倍體樣本 手里有參考基因組你就具備了用 Snippy 的全部前提。? 問題五怎么用倉庫自帶測試數(shù)據(jù)跑出第一份 SNP 檢出倉庫的test目錄里躺著三件套example.fna參考序列、example.gbk帶注釋的參考、example.bed區(qū)域文件。就拿它們當(dāng)?shù)谝环菥毷植牧稀U鎸?shí) reads 文件太大先用 wgsim 從參考序列模擬一對雙端 reads這也是官方測試流程的慣用做法wgsim -S 1 -h -r 0.005 -N 12000 -1 100 -2 100 -d 200 example.fna R1.fq R2.fq上面參數(shù)的意思隨機(jī)變異率 0.5%生成 12000 對長度 100 bp 的 reads插入片段約 200 bp。然后正式開跑snippy --cpus 4 --outdir my_first_run --ref example.fna --R1 R1.fq --R2 R2.fq日志里會依次出現(xiàn) bwa 比對、freebayes 變異識別等階段最后以類似這樣的三行收尾Walltime used: 3 min, 42 sec Results folder: my_first_run Done.進(jìn)目錄看看都產(chǎn)出了什么ls my_first_run你會拿到一整套文件snps.vcf標(biāo)準(zhǔn) VCF 變異文件、snps.tab易讀表格、snps.bam比對文件、snps.gff、snps.bed、snps.html還有snps.consensus.fa——把變異全部回貼到參考序列上得到的修正版基因組。用head -5 my_first_run/snps.tab看表格前幾行先認(rèn)識六列核心信息CHROM變異所在的序列名POS位置TYPE變異類型常見 snp / mnp / ins / del / complexREF參考堿基ALT樣本里由 reads 支持的堿基EVIDENCE支持各堿基的 reads 計數(shù)如果你把--ref換成example.gbk這種帶注釋的文件表格還會多出一串基因相關(guān)列FTYPE、STRAND、GENE、PRODUCT以及 snpEff 預(yù)測的EFFECT。變異落在哪個基因、是不是錯義突變、會不會改變氨基酸一眼就能看到——這是 Snippy 很貼心的設(shè)計。帶走什么到這你手里已經(jīng)有第一張變異表格了后面所有進(jìn)階玩法都建立在這套產(chǎn)物上。 問題六十幾個樣本要跟同一參考比對怎么批量做樣本少時一個個手敲snippy還湊合一多就是純折磨。Snippy 為此準(zhǔn)備了批量入口snippy-multi你只要準(zhǔn)備一個制表符分隔的清單文件input.tab一行一個樣本Isolate1 /path/to/R1.fq.gz /path/to/R2.fq.gz Isolate2 /path/to/SE.fq.gz Isolate3 /path/to/contigs.fa同一行內(nèi)用 Tab 分隔雙端、單端、contigs 三種形態(tài)可以混排。然后生成并檢查批量腳本snippy-multi input.tab --ref Reference.gbk --cpus 16 runme.sh less runme.sh確認(rèn)腳本內(nèi)容沒問題再放行sh runme.sh跑完后每個樣本有獨(dú)立結(jié)果目錄而且snippy-multi會在最后自動調(diào)用snippy-core把每個樣本都有覆蓋的基因組位置挑出來做核心基因組比對產(chǎn)出core.aln多序列比對、core.vcf多樣本 VCF等文件。收尾時你會看到類似這樣的匯總Found 2814 core SNPs from 96615 SNPs.意思是全部 96615 個變異位點(diǎn)里有 2814 個是各樣本都覆蓋到的核心 SNP。把core.aln丟給 FastTree 這類建樹工具就能畫系統(tǒng)發(fā)育樹。帶走什么一個清單文件加一條命令Snippy 批量變異檢測和核心基因組比對就都齊了。 問題七真實(shí)數(shù)據(jù)不按套路出牌四個高頻麻煩怎么解測試數(shù)據(jù)很乖真實(shí)數(shù)據(jù)可不一定。下面四個場景按出現(xiàn)頻率排遇到直接對照處理。場景 1報command not found多半是 PATH 沒配好。用which snippy、which bwa逐個排查哪個找不到就把它所在目錄加進(jìn) PATH或者干脆改用絕對路徑。場景 2深度太高跑得特別慢有的樣本測了幾百上千倍深度而大部分變異在 50~100x 深度就能可靠檢出。此時按比例隨機(jī)抽讀即可比如 1000x 想降到 100xsnippy --subsample 0.1 --outdir out --ref ref.fna --R1 R1.fq.gz --R2 R2.fq.gz日志里出現(xiàn) Sub-sampling reads at rate 0.1 就說明生效了提速立竿見影。場景 3只想篩特定區(qū)域的突變比如只關(guān)心耐藥基因上的變異那就把感興趣的區(qū)域?qū)戇M(jìn) BED 文件用--targets限定范圍能省下大量計算時間snippy --targets sites.bed --outdir out --ref ref.fna --R1 R1.fq.gz --R2 R2.fq.gz場景 4手里只有 contigs原始 reads 早就丟了別慌把 contigs 文件直接交給--ctgsSnippy 會把它拆成 250 bp 的合成短讀再做比對snippy --outdir mut1 --ref ref.gbk --ctgs mut1.fasta它產(chǎn)出的結(jié)果目錄和 reads 樣本完全兼容可以混在一起參與snippy-core分析。帶走什么這四個場景覆蓋了新手階段九成以上的跑不動逐個對照就能解。 收尾裝好之后順手做掉這三件小事最后給你三個動作花不了多少時間卻能讓后面的路順暢很多把測試數(shù)據(jù)完整走一遍snippy --check → 單樣本 → 批量全流程把每步的預(yù)期輸出記在腦子里再碰真實(shí)數(shù)據(jù)。真實(shí)樣本開跑前備份原始 reads給每個樣本起一個清晰可追溯的 ID。記下版本號。變異檢測結(jié)果跟版本強(qiáng)相關(guān)寫文章或匯報時附上snippy --version的輸出結(jié)果才可復(fù)現(xiàn)、可信。從裝環(huán)境到批量跑完 SNP 檢出這條鏈路已經(jīng)全部打通。接下來無論是十幾個樣本的群體分析還是基于核心 SNP 比對建系統(tǒng)發(fā)育樹你手里的數(shù)據(jù)地基都已經(jīng)足夠牢靠。【免費(fèi)下載鏈接】snippy:scissors: :zap: Rapid haploid variant calling and core genome alignment項目地址: https://gitcode.com/gh_mirrors/sn/snippy創(chuàng)作聲明:本文部分內(nèi)容由AI輔助生成(AIGC),僅供參考