做單細(xì)胞的人繞不開兩件事一是在 NCBI 上把別人發(fā)表過的 10X 單細(xì)胞測序原始數(shù)據(jù)扒下來二是把它喂進(jìn) Cell Ranger 跑出表達(dá)矩陣。論文里那句輕描淡寫的raw data are available in GEO under accession GSEXXXXXX背后其實(shí)是一整套從 SRA 到 FASTQ、從文件命名到樣本表對齊的體力活。我前前后后處理過三十多個(gè)公開數(shù)據(jù)集從人的 PBMC 到小鼠的腦組織踩過的坑包括 barcode 被當(dāng)成普通 read 丟掉、文件名讓 Cell Ranger 直接報(bào) sample 找不到、下到一半磁盤爆掉、跑了兩天的結(jié)果質(zhì)控一塌糊涂。這篇就把 NCBI 下載 10X 單細(xì)胞測序原始數(shù)據(jù)到 Cell Ranger 分析這條完整鏈路拆開來聊從檢索、下載、格式轉(zhuǎn)換、輸入準(zhǔn)備一直到運(yùn)行和排錯(cuò)盡量把每一步的為什么也講清楚。不管你是剛進(jìn)實(shí)驗(yàn)室的研一新生還是已經(jīng)會(huì)跑 Cell Ranger 但每次下數(shù)據(jù)都手忙腳亂的老手這套流程應(yīng)該都能直接用。1. 這條鏈路到底長什么樣先拆解再動(dòng)手很多人第一次從 NCBI 下 10X 數(shù)據(jù)是直接在 GEO 頁面找Download按鈕結(jié)果點(diǎn)進(jìn)去發(fā)現(xiàn)全是 .sra、.bam、.tar 之類的文件根本不知道怎么變成 Cell Ranger 要的 fastq 目錄。問題的根源在于GEO、SRA 和 Cell Ranger 三者對原始數(shù)據(jù)的理解完全不一樣。1.1 為什么 NCBI 下載這件事比想象中麻煩GEO 是一個(gè)數(shù)據(jù)集展示層它關(guān)心的是樣本分組、實(shí)驗(yàn)設(shè)計(jì)、表型信息SRA 是一個(gè)原始序列倉庫它只負(fù)責(zé)把測序儀吐出來的堿基流保存成一種緊湊的二進(jìn)制格式而 Cell Ranger 只認(rèn)一種東西命名符合規(guī)范、read 結(jié)構(gòu)正確的 gzip 壓縮 fastq。這三個(gè)體系之間的轉(zhuǎn)換沒有任何一個(gè)按鈕能一鍵完成。更麻煩的是 10X 數(shù)據(jù)本身的特殊性。常規(guī)的 WGS 或 RNA-seq 是雙端測序R1 和 R2 各自就是一條真實(shí)的序列片段但 10X Chromium 的 read 結(jié)構(gòu)是錯(cuò)位的——R1 裝的是 16bp barcode 加 10bp UMIv2 版本 26bpv3 版本 28bpR2 才是真正的 cDNA 序列中間還可能夾著一條 I1 的 sample index。如果下載工具按常規(guī)的雙端邏輯處理極容易把 barcode 那條 read 直接扔掉。我見過最典型的翻車場景是這樣的用fastq-dump不帶任何參數(shù)跑出來兩個(gè)文件扔進(jìn) Cell Ranger 報(bào)錯(cuò)Read 1 is too short因?yàn)橄螺d工具默認(rèn)只取了 biological reads把 barcode 所在的 index read 當(dāng)成了技術(shù)性 read 過濾掉了。搞清楚這一點(diǎn)后面所有操作就都有了依據(jù)。1.2 三種常見的數(shù)據(jù)交付形態(tài)在 GEO 上翻數(shù)據(jù)集的時(shí)候原始數(shù)據(jù)一般以三種形態(tài)存在處理難度差別很大形態(tài)典型文件處理難度推薦做法交付 processed 矩陣.h5、.mtx、.tsv最低直接用 Seurat/Scanpy 讀不必跑 Cell Ranger交付 bam 文件possorted_genome_bam.bam中等用cellranger bamtofastq反推 fastq只掛 SRA 號(hào)SRR/SRP/PRJNA最高SRA Toolkit 或 ENA 下載第一種情況最省事很多 2019 年之后發(fā)表的文章會(huì)把 filtered_feature_bc_matrix.h5 直接放在 Supplementary Files 里。但注意這類矩陣往往已經(jīng)被作者做過一輪過濾你想復(fù)現(xiàn)原文的 Cell Ranger 參數(shù)或者重新做 QC 就無能為力了。真正要做方法學(xué)復(fù)現(xiàn)、要自己調(diào)--expect-cells這類參數(shù)的還是得拿原始 fastq。1.3 完整鏈路與耗時(shí)預(yù)估把整條鏈路鋪開大致是這么幾步在 GEO 定位到研究數(shù)據(jù)集順著鏈接找到對應(yīng)的 SRA/BioProject accession用 SRA Toolkit 或者 ENA 下載成 fastq檢查 read 長度確認(rèn) barcode 在不在按 Cell Ranger 規(guī)則重命名準(zhǔn)備參考基因組最后cellranger count跑完拿矩陣。時(shí)間上一個(gè)常規(guī)的 10X 3 樣本約 4 億 reads下載取決于網(wǎng)絡(luò)和源站ENA 走直連大概半小時(shí)到兩小時(shí)fasterq-dump 轉(zhuǎn)換比下載更慢磁盤 IO 是瓶頸20 到 60 分鐘Cell Ranger count 在 16 核 64G 的機(jī)器上大約 2 到 5 小時(shí)。也就是說單個(gè)樣本從零到矩陣留出 8 小時(shí)的窗口比較穩(wěn)妥。批量處理 10 個(gè)樣本的話建議頭一天晚上掛上下載第二天白天跑轉(zhuǎn)換和比對排期才不至于太緊張。2. 在 NCBI 上精準(zhǔn)定位可用的 10X 數(shù)據(jù)集下載慢一點(diǎn)還能忍最虧的是辛辛苦苦下完才發(fā)現(xiàn)數(shù)據(jù)集根本不適合自己的問題。檢索階段多花二十分鐘能省掉后面一整天的返工。2.1 GEO 檢索從關(guān)鍵詞到 accession 號(hào)打開 GEO DataSets 的檢索框不要只輸入研究主題要把平臺(tái)特征一起加上。我常用的組合是(10x OR Chromium OR single cell RNA) AND 你的組織或疾病關(guān)鍵詞。括號(hào)和 OR 的寫法在 GEO 里是支持的能比空格分詞精準(zhǔn)不少。拿到候選 GSE 號(hào)之后進(jìn)詳情頁重點(diǎn)看三處頂部的 Platform 字段、Samples 列表里的 Organism 和 Source Name、底部的 Series 關(guān)聯(lián)的 BioProject。Platform 如果是 GPL24676 或者 GPL18573 這類帶 10x 字樣的基本可以確定是 Chromium 3 或 5 試劑盒如果 Platform 描述里出現(xiàn)10x Genomics Visium那是空間轉(zhuǎn)錄組跟單細(xì)胞完全是兩套流程別下錯(cuò)了。從 GSE 頁面左下角的Relations里點(diǎn)進(jìn) BioProject再到 SRA Run Selector 頁面才能看到真正的測序數(shù)據(jù)。SRA Run Selector 是我最推薦的中轉(zhuǎn)站它把每個(gè) run 的樣本名、reads 數(shù)、堿基量、文件大小都列成一張可下載的表格還能直接勾選導(dǎo)出 Accession List 和 Metadata 兩份 txt。2.2 判斷數(shù)據(jù)集是不是 Chromium 平臺(tái)這里有個(gè)血淚教訓(xùn)GEO 上標(biāo)的single cell不一定就是 10X早期用 Smart-seq2、CEL-Seq、Drop-seq 平臺(tái)做的也是單細(xì)胞這些數(shù)據(jù)的 read 結(jié)構(gòu)和 10X 完全不同喂給 Cell Ranger 會(huì)直接崩。判斷方法有三條看 README 或者 Methods 部分有沒有出現(xiàn)Chromium Controller、10x Genomics字樣看 SRA 里 run 的 read 數(shù)量10X 一個(gè) run 動(dòng)輒上億Smart-seq2 單細(xì)胞通常只有幾百萬看 read 長度10X 的 R1 通常是 26bp 或 28bp這個(gè)長度非常具有辨識(shí)度。如果三條里兩條對不上先別急著下用fastq-dump -X 10 --split-files抓前 10 條序列看一眼長度幾秒鐘的事比下完 40G 再發(fā)現(xiàn)不對劃算太多。2.3 樣本數(shù)與元信息核對清單決定開下之前我會(huì)在表格里核對這幾項(xiàng)run 總數(shù)和樣本數(shù)是否一致有些數(shù)據(jù)集一個(gè)樣本拆成了多個(gè) lane需要合并每個(gè) run 的 spots 數(shù)也就是細(xì)胞數(shù)有沒有明顯偏小低于 500 個(gè)細(xì)胞的樣本后續(xù)分析價(jià)值有限物種和參考基因組版本人、鼠、還是混合物種是否用了 GRCh38 還是 hg19組織處理方式是新鮮樣本還是凍存是細(xì)胞懸液還是細(xì)胞核后者在做 QC 時(shí)閾值要放寬是否有配套的 Cell Ranger 輸出文件有的話可以用來做下游結(jié)果的交叉驗(yàn)證。把這些信息整理成一張表后面構(gòu)建樣本表和排查問題時(shí)隨時(shí)能翻比反復(fù)回 GEO 頁面找要快得多。3. 下載環(huán)節(jié)SRA Toolkit 與 ENA 兩條路怎么選真正開始下載其實(shí)有兩條成熟路線一是用 SRA Toolkit 從 NCBI 官方源拉二是繞道 ENA 直接下 fastq。兩條路我都長期在用各有適用場景。3.1 SRA Toolkit 安裝與初始配置首推 conda 安裝避免自己編譯的麻煩conda create -n sra python3.10 conda activate sra conda install -c bioconda sra-tools pigz裝完之后立刻做一件事把默認(rèn)的緩存目錄改掉。SRA Toolkit 默認(rèn)把數(shù)據(jù)放在~/ncbi/public/sra/這個(gè)目錄直接吃 home 分區(qū)空間家里分區(qū)一滿整個(gè)系統(tǒng)都可能出問題。在~/.ncbi/user-settings.mkfg里加兩行echo /repository/user/main/public/root /data/sra_cache ~/.ncbi/user-settings.mkfg echo /repository/user/main/public/cache-enabled true ~/.ncbi/user-settings.mkfg另外建議把prefetch的最大文件限制調(diào)大否則遇到大樣本會(huì)直接拒絕下載vdb-config --set /sratoolkit/refseq/download-max-size 100000000000這些配置看著瑣碎但每一條都是我實(shí)際踩坑之后補(bǔ)上的——緩存目錄設(shè)在 home 分區(qū)被撐爆、大文件被默認(rèn)限制卡住都是新手最常見的兩種翻車方式。3.2 prefetch 與 fasterq-dump 的標(biāo)準(zhǔn)組合正式下載分兩步走先用prefetch把 .sra 二進(jìn)制文件拉到本地再用fasterq-dump轉(zhuǎn)換成 fastq。之所以不一步到位是因?yàn)?prefetch 支持?jǐn)帱c(diǎn)續(xù)傳fasterq-dump 不支持網(wǎng)絡(luò)一抖動(dòng)前功盡棄。# 第一步下載 sra 文件支持?jǐn)帱c(diǎn)續(xù)傳 prefetch SRR1234567 -O /data/sra_cache批量下載的話可以把 accession 寫進(jìn) txt 每個(gè)一行然后cat srr_list.txt | while read id; do prefetch $id -O /data/sra_cache; done第二步轉(zhuǎn)換關(guān)鍵是加對參數(shù)fasterq-dump /data/sra_cache/SRR1234567.sra \ --split-files \ --include-technical \ --threads 8 \ --temp /data/tmp \ -O /data/fastq/SRR1234567--split-files保證 paired reads 被拆成獨(dú)立的文件這是 Cell Ranger 的硬性要求--include-technical是 10X 數(shù)據(jù)的救命參數(shù)加了它才會(huì)把 index read 一起輸出。如果這個(gè) run 確實(shí)沒有保存 index加了也不會(huì)有副作用只是不會(huì)多出第三個(gè)文件而已。--temp指定臨時(shí)目錄也很重要fasterq-dump 生成過程中臨時(shí)文件體積能到最終結(jié)果的近兩倍默認(rèn)走 /tmp 很容易把系統(tǒng)盤塞滿。轉(zhuǎn)換完的 fastq 默認(rèn)不壓縮逐個(gè)用 pigz 壓一下pigz -p 8 /data/fastq/SRR1234567/*.fastqpigz 是多線程版的 gzip同樣大小的文件壓縮速度比單線程 gzip 快好幾倍這個(gè)細(xì)節(jié)在大批量處理時(shí)省下的時(shí)間相當(dāng)可觀。3.3 ENA 直連下載更省事的一條捷徑如果只是想拿到 fastq其實(shí) ENA歐洲核酸檔案庫比 NCBI 貼心很多它把 SRA 里的數(shù)據(jù)預(yù)先轉(zhuǎn)成了 fastq.gz還拆好了文件、起好了名字。拿到 BioProject 號(hào)之后先拉一份文件清單curl https://www.ebi.ac.uk/ena/portal/api/filereport?accessionPRJNA123456resultread_runfieldsrun_accession,fastq_ftp,fastq_bytes,read_countformattsv ena_manifest.tsv打開這個(gè) tsvfastq_ftp一列就是下載地址。批量下載很簡單tail -n 2 ena_manifest.tsv | cut -f2 | tr ; \n | sed s|^|https://| urls.txt aria2c -i urls.txt -j 4 -x 8 -d /data/fastq/enaENA 的優(yōu)勢非常明顯免去了 fasterq-dump 那一步源站帶寬也普遍比 NCBI 更友好實(shí)測同樣的數(shù)據(jù)集速度能差三到五倍。但它也不是萬能——有些提交者上傳的 SRA 數(shù)據(jù)經(jīng)過了額外處理ENA 上未必 100% 還原原始 read 結(jié)構(gòu)所以下完仍然要做第 4 章的長度檢查。我的習(xí)慣是ENA 有就直接用 ENAENA 缺數(shù)據(jù)再回 SRA Toolkit。4. 整理成 Cell Ranger 能吃的輸入下載完成只是拿到了原材料真正決定成敗的是把它們整理成 Cell Ranger 認(rèn)識(shí)的樣子。這一步不出錯(cuò)后面的比對基本就穩(wěn)了。4.1 先看清 read 結(jié)構(gòu)barcode 到底在哪個(gè)文件里拿到一堆 fastq 之后第一件事是抽前幾條序列看長度head -n 8 /data/fastq/SRR1234567/SRR1234567_1.fastq head -n 8 /data/fastq/SRR1234567/SRR1234567_2.fastq對照下表判斷平臺(tái)版本R1 長度R2 長度是否有 I1Chromium 3 v22698可選Chromium 3 v32891 或 150可選Chromium 5 v1/v22698可選Chromium 3 v42890可選如果 _1 文件里長度是 26 或 28說明 barcode 結(jié)構(gòu)完好如果 _1 的長度接近 100 甚至更長那大概率是提交者把 barcode 拼到了 cDNA 前面或者數(shù)據(jù)被重新處理過這時(shí)候需要跟原始發(fā)表文章核實(shí)必要時(shí)聯(lián)系作者要原始數(shù)據(jù)。如果發(fā)現(xiàn)根本沒出 _1 文件只有 _2那說明--include-technical沒加對或者這個(gè) SRA 條目本身沒保留 index。還有一個(gè)隱形的坑有些數(shù)據(jù)集的 _1 長度正確但用錯(cuò)了 read 方向測序時(shí)把 cDNA 放在了 _1 里。這種情況很少但存在判斷方法是看 _1 的堿基分布barcode 那部分因?yàn)橐淙?whitelist四種堿基的分布會(huì)比較均勻而 cDNA 有明顯偏向性。4.2 文件重命名規(guī)則與 sample 參數(shù)對齊Cell Ranger 對文件名的要求非常死板必須符合這個(gè)模式[SAMPLE_NAME]_S[NUM]_L[LANE]_R[1-2]_[NUM].fastq.gz從 SRA 出來的是SRR1234567_1.fastq.gz和SRR1234567_2.fastq.gz完全不匹配。重命名示例cd /data/fastq/SRR1234567 mv SRR1234567_1.fastq.gz Pbmc_Donor1_S1_L001_R1_001.fastq.gz mv SRR1234567_2.fastq.gz Pbmc_Donor1_S1_L001_R2_001.fastq.gz這里的Pbmc_Donor1就是后面--sample參數(shù)要填的值必須完全一致。我強(qiáng)烈建議用見名知意的樣本名而不是 SRR 號(hào)因?yàn)楹竺娉鰣D的時(shí)候軸標(biāo)簽直接就是這個(gè)名字用 SRR 號(hào)看著非常難受。如果同一個(gè)樣本被拆成多個(gè) run比如 SRR111111、SRR222222 屬于同一個(gè)樣本重命名時(shí)把它們的 sample 部分寫成一致S 和 L 編號(hào)不同即可Pbmc_Donor1_S1_L001_R1_001.fastq.gz Pbmc_Donor1_S1_L002_R1_001.fastq.gzCell Ranger 會(huì)自動(dòng)識(shí)別屬于同一個(gè) sample 的多個(gè) lane 并合并處理不用手動(dòng) cat。4.3 參考基因組下載與版本選擇Cell Ranger 需要 10x 官方預(yù)構(gòu)建的參考基因組不是隨便一個(gè) FASTA 加 GTF。下載地址在 10x Genomics 支持頁面的Reference Genomes部分常見的兩類構(gòu)建refdata-gex-GRCh38-2020-A人類2020 版最通用refdata-gex-mm10-2020-A小鼠對應(yīng) GRCm38/mm10。下載的是一個(gè)約 11G 的 tar.gzwget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz tar -xzvf refdata-gex-GRCh38-2020-A.tar.gz -C /data/refs/版本選擇上2020-A 是經(jīng)過最多驗(yàn)證的版本除非原始文章明確說自己用了舊版否則優(yōu)先選它。如果用 GRCh38-2020-A 而原文用的是 hg19 或 GRCh38-1.2.0基因注釋和基因名會(huì)有細(xì)微差別做出來的下游結(jié)果和原文對不上號(hào)不能算 bug只能說是版本差異這一點(diǎn)寫方法學(xué)的時(shí)候要說清楚。5. Cell Ranger count 實(shí)操與資源規(guī)劃輸入都備好之后就進(jìn)入真正跑分析的環(huán)節(jié)。這個(gè)環(huán)節(jié)最貴的是時(shí)間參數(shù)設(shè)錯(cuò)一次可能白跑半天所以值得把每個(gè)參數(shù)都過一遍。5.1 命令逐參數(shù)拆解cellranger count \ --idPbmc_Donor1_run1 \ --transcriptome/data/refs/refdata-gex-GRCh38-2020-A \ --fastqs/data/fastq \ --samplePbmc_Donor1 \ --expect-cells5000 \ --localcores16 \ --localmem64 \ --create-bamtrue每個(gè)參數(shù)的作用和踩坑點(diǎn)--id輸出目錄名如果目錄已存在會(huì)直接報(bào)錯(cuò)想重跑要加--force或者換個(gè) id別傻乎乎刪了又重跑。--fastqsfastq 所在目錄Cell Ranger 會(huì)自動(dòng)遞歸查找所以目錄層級(jí)不必完全精確但樣本名匹配必須精確。--sample這是最容易出錯(cuò)的地方它必須和文件名_S1_前面的部分完全一致。大小寫、下劃線數(shù)量都要對上差一個(gè)字符就會(huì)報(bào)no input FASTQs were found。--expect-cells告訴算法預(yù)期細(xì)胞數(shù)。10X 官方的自動(dòng)估計(jì)在細(xì)胞數(shù)差異大的樣本上不太靠譜如果已知樣本上機(jī)時(shí)目標(biāo)細(xì)胞數(shù)直接填比如說 5000 或 10000如果不知道可以先不填跑一次看估計(jì)值再補(bǔ)跑。--localcores和--localmem必須顯式寫否則 Cell Ranger 會(huì)按機(jī)器全部核數(shù)去要資源在共享服務(wù)器上很容易把別人擠下去也容易被系統(tǒng) OOM killer 干掉。--create-bam默認(rèn) true輸出 bam 文件占空間很大如果只想要矩陣可以設(shè)成 false 省空間。5.2 時(shí)間與內(nèi)存的實(shí)際估算Cell Ranger count 是內(nèi)存密集型內(nèi)存需求的經(jīng)驗(yàn)公式大致是每 1000 個(gè)細(xì)胞配 1G 內(nèi)存但對于 reads 數(shù)非常高的樣本這個(gè)公式會(huì)低估。實(shí)際經(jīng)驗(yàn)值細(xì)胞數(shù)reads 數(shù)推薦內(nèi)存16 核耗時(shí)30001 億32G1.5-2 小時(shí)50002 億48G2-3 小時(shí)100004 億64G3-5 小時(shí)200008 億96-128G6-10 小時(shí)內(nèi)存不夠的典型癥狀是任務(wù)在中途毫無征兆地被 kill日志里只有一行 Killed。這種情況下先看系統(tǒng) OOM 日志然后按上表加內(nèi)存再跑。共享服務(wù)器上跑之前最好先free -h看一眼可用內(nèi)存以及用squeue或htop確認(rèn)沒人跟你搶。5.3 輸出目錄與質(zhì)控指標(biāo)怎么看跑完之后輸出目錄里重要的東西有這些outs/web_summary.html可視化報(bào)告第一個(gè)該看的東西outs/metrics_summary.csv關(guān)鍵質(zhì)控指標(biāo)的表格版outs/filtered_feature_bc_matrix/過濾后的表達(dá)矩陣下游分析入口outs/raw_feature_bc_matrix/原始矩陣包含空液滴做 SoupX 之類的背景校正會(huì)用到outs/possorted_genome_bam.bam排序后的比對結(jié)果占空間最多。web_summary 里最該盯的幾個(gè)數(shù)字和大致閾值指標(biāo)理想值說明Estimated Number of Cells與預(yù)期接近差一個(gè)數(shù)量級(jí)要警惕Median Genes per Cell 500低于 200 一般說明樣本質(zhì)量差Median UMI per Cell 1000低于 500 得上游找原因Reads Mapped Confidently to Transcriptome 60%低于 50% 檢查參考版本Fraction Reads in Cells 70%偏低說明空液滴占比高Q30 Bases in Barcode 85%偏低是測序問題得換數(shù)據(jù)實(shí)際評估的時(shí)候我不會(huì)只看單項(xiàng)而是把幾個(gè)指標(biāo)一起看。比如同時(shí)出現(xiàn)細(xì)胞數(shù)正常但 median genes 只有 300和Fraction Reads in Cells 只有 50%那大概率是樣本本身 RNA 降解嚴(yán)重得回原文確認(rèn)他們用的是什么處理方法是不是冷凍組織或者細(xì)胞核樣本后者的中位基因數(shù)本來就會(huì)偏低一些。6. 報(bào)錯(cuò)排查與實(shí)操心得跑流程的過程中報(bào)錯(cuò)是常態(tài)。下面把遇到的典型問題和處理方式整理出來遇到問題直接查表能省不少時(shí)間。6.1 下載與轉(zhuǎn)換階段的常見問題問題一prefetch 報(bào)maximum file size exceeded。原因是 SRA Toolkit 默認(rèn)的文件大小限制太小解決方案是在配置里把download-max-size調(diào)大具體命令前面 3.1 節(jié)已經(jīng)給了。問題二fasterq-dump 跑到一半提示磁盤空間不足。fasterq-dump 在轉(zhuǎn)換過程中會(huì)產(chǎn)生幾乎和最終結(jié)果等體積的臨時(shí)文件然后把它們重新拼裝。實(shí)際預(yù)留的空間要比最終結(jié)果大兩倍以上。--temp參數(shù)一定要顯式指定到空間充足的分區(qū)。問題三轉(zhuǎn)換完發(fā)現(xiàn)只有一個(gè) fastq 文件。說明這個(gè) run 可能本身就是單端或者--split-files沒加或者數(shù)據(jù)在 SRA 里的元信息把兩條 read 標(biāo)記成了同一條 spot。先加--split-spot試試實(shí)在不行用fastq-dump老命令加--split-files再跑一次兩個(gè)工具的內(nèi)部邏輯不完全一樣。問題四ENA 下載下來的 fastq.gz 解壓后發(fā)現(xiàn)文件截?cái)?。多半是下載中斷了aria2c支持?jǐn)帱c(diǎn)續(xù)傳重跑一次同樣的命令會(huì)接著下。如果重跑后還是壞的用gzip -t校驗(yàn)一遍確認(rèn)完整性再往下一步走。6.2 運(yùn)行階段的常見問題問題一The chemistry of this run was not detected。Cell Ranger 無法從 read 長度自動(dòng)推斷試劑盒版本。這通常發(fā)生在 read 長度不太標(biāo)準(zhǔn)的情況下。解決辦法是手動(dòng)指定--chemistrySC3Pv3之類的參數(shù)具體是 v2 還是 v3 靠 R1 是 26 還是 28 判斷。問題二no input FASTQs were found for sample xxx。回到第 4.2 節(jié)檢查三件事文件后綴是不是 .fastq.gz不能是 .fq.gz、文件名里的_R1__R2_有沒有漏、--sample參數(shù)和文件名前綴是否完全一致。我遇到過一次是文件名里用了中文下劃線肉眼看不出來折騰了一小時(shí)才發(fā)現(xiàn)。問題三Cell Ranger 跑到比對階段內(nèi)存被 kill。如果是集群環(huán)境先看 SLURM 或 PBS 的資源申請是否設(shè)置對如果是本地機(jī)器看dmesg | grep -i killed確認(rèn)是不是 OOM。加內(nèi)存或者先用--subsample參數(shù)小規(guī)模跑一遍驗(yàn)證流程。問題四web_summary 里 Estimated Number of Cells 只有幾百但實(shí)際材料明明有幾千個(gè)細(xì)胞。兩個(gè)常見原因一是--expect-cells沒有設(shè)置或者默認(rèn)值太低導(dǎo)致算法偏向保守二是原始數(shù)據(jù)質(zhì)量問題RNA 含量低導(dǎo)致細(xì)胞識(shí)別不出來。可以設(shè)一個(gè)明確的--expect-cells重跑一次看看結(jié)果是否變化。6.3 幾張私藏的經(jīng)驗(yàn)卡片用了這么多年攢下這么幾條不太寫在文檔里的心得都用得上?？ㄆ挥肋h(yuǎn)保留一套原始 SRA 文件。fasterq-dump 轉(zhuǎn)換一次可能要幾小時(shí)如果后面發(fā)現(xiàn)參數(shù)錯(cuò)了要重新轉(zhuǎn)換沒有原始 .sra 就得重新下載。我習(xí)慣在 /data/sra_cache 里留一份 .sra等所有分析都確認(rèn)沒問題了再統(tǒng)一清理。卡片二樣本名里不要出現(xiàn)空格和特殊字符。Cell Ranger 對文件名里的特殊字符非常敏感用下劃線連接是最安全的。樣本名里的群體、處理批次信息用Ctrl_、Treat_這種前綴別用C-1、T#2這種符號(hào)。卡片三md5 校驗(yàn)一定要做。大文件下載完之后跑一次md5sum -c花不了幾分鐘但能避免后面在錯(cuò)誤數(shù)據(jù)上浪費(fèi)一整天。尤其是走 ENA 或者第三方鏡像的時(shí)候文件缺失或者截?cái)嗟母怕什⒉坏?。卡片四下游分析之前先做一次聚類目視檢查。拿到 filtered_feature_bc_matrix 之后別急著做差異分析先跑個(gè) PCA 和 UMAP 看一眼有沒有明顯異常樣本。有些時(shí)候 Cell Ranger 的指標(biāo)看起來正常但實(shí)際細(xì)胞群的分布非常奇怪這時(shí)候就該回上游重新審視參數(shù)設(shè)置。卡片五方法學(xué)描述要寫清楚版本號(hào)。NCBI、SRA Toolkit、Cell Ranger、參考基因組這四樣?xùn)|西的版本號(hào)都要記下來投稿的時(shí)候?qū)徃迦撕芸赡軙?huì)問。Cell Ranger 在 6.0 之后對 UMI 的處理有一處調(diào)整不同版本跑出來的矩陣在個(gè)別基因上會(huì)有差異寫清楚版本號(hào)既是對讀者負(fù)責(zé)也是避免日后自己都說不清楚。最后再補(bǔ)一個(gè)我用了很多年的小技巧把整個(gè)流程寫成一個(gè) shell 腳本把 accession 號(hào)做成參數(shù)傳進(jìn)去從 prefetch 到 cellranger count 一條龍跑完。第一個(gè)樣本手工走通之后后面幾十個(gè)樣本照著腳本批量跑犯錯(cuò)的概率會(huì)大幅下降也方便日后復(fù)現(xiàn)。腳本里每跑完一步就寫一行日志到文件出了問題回頭翻日志定位比盯著終端輸出回溯要高效得多。