企業(yè)營銷網(wǎng)站建設(shè)費用,商務(wù)網(wǎng)站管理與建設(shè),深圳龍崗現(xiàn)在算什么風(fēng)險地區(qū),wordpress返回首頁一直加載全基因組重測序上游分析流程#xff5c;從軟件部署到變異檢測#xff0c;超細致實操指南 作為科研新手#xff0c;第一次上手全基因組重測序數(shù)據(jù)處理時#xff0c;我踩過不少軟件安裝的坑、碰過參數(shù)設(shè)置的雷。如今整理出這份超詳細流程#xff0c;從前期準備到最終變異過…全基因組重測序上游分析流程從軟件部署到變異檢測超細致實操指南作為科研新手第一次上手全基因組重測序數(shù)據(jù)處理時我踩過不少軟件安裝的坑、碰過參數(shù)設(shè)置的雷。如今整理出這份超詳細流程從前期準備到最終變異過濾每一步都標注了關(guān)鍵注意事項跟著練一遍就能快速上手。覺得有用的話別忘了點贊收藏哦適用場景動植物全基因組重測序上游分析變異檢測核心流程核心工具BWA、Samtools、Picard、GATK4全開源附conda安裝命令0 前期準備軟件部署與環(huán)境搭建重測序上游分析的核心工具就4個全部開源可通過conda快速安裝推薦Linux集群或Mac OS系統(tǒng)Windows建議用WSL2。重點注意GATK4的版本特性——雖然是新版本但100%開源且適配大規(guī)模數(shù)據(jù)是未來主流本文全程基于GATK4構(gòu)建流程。工具名稱核心功能conda安裝命令注意事項BWANGS數(shù)據(jù)與參考基因組比對conda install -c bioconda bwaC語言編寫需系統(tǒng)支持編譯SamtoolsBAM/SAM文件處理排序、索引等conda install -c bioconda samtools與Picard功能互補必裝工具Picard標記重復(fù)序列、文件格式處理conda install -c bioconda picardJava編寫需Java 1.8環(huán)境GATK4變異檢測、基因型推斷conda install -c bioconda gatk44.x比3.x更適配集群全開源避坑指南安裝生信軟件務(wù)必加-c bioconda指定頻道避免下載到舊版本或錯誤包。Java版本過低會導(dǎo)致Picard報錯可通過java -version檢查低于1.8則需重新安裝。1 原始數(shù)據(jù)質(zhì)控快速過一遍心里有底現(xiàn)在測序公司交付的基本都是經(jīng)過初步處理的clean data但自己驗證一步更放心。核心用兩個工具FastQC可視化展示數(shù)據(jù)質(zhì)量堿基分布、測序錯誤率等命令fastqc read_1.fq.gz read_2.fq.gzfastp批量清洗數(shù)據(jù)去除接頭、低質(zhì)量堿基命令fastp -i read_1.fq.gz -I read_2.fq.gz -o clean_1.fq.gz -O clean_2.fq.gz如果FastQC報告中“Per base sequence quality”出現(xiàn)紅色區(qū)域或接頭污染率高就需要用fastp調(diào)整參數(shù)如--cut_front去除前端低質(zhì)量堿基再處理。2 核心流程從序列比對到變異檢測這部分是重測序分析的核心每一步都有明確的邏輯和避坑點跟著步驟走準沒錯。2.1 序列比對給短讀長“找家”NGS測出來的短序列read是隨機打亂的必須通過比對找到它們在參考基因組上的位置。BWA是目前最權(quán)威的工具核心靠“索引構(gòu)建比對”兩步。步驟1構(gòu)建參考基因組索引索引能讓BWA快速定位序列相當(dāng)于給參考基因組建“目錄”。命令超簡單bwa index genome.fasta運行后會生成5個以genome.fasta為前綴的文件.amb/.ann/.bwt等這些是比對的關(guān)鍵別刪步驟2雙末端序列比對重測序常用雙末端測序PE兩個fq文件分別對應(yīng)DNA片段的兩端比對時要一起輸入。這里有個超級關(guān)鍵的參數(shù)——Read Group-R直接影響后續(xù)GATK分析bwa mem -t 4 -R RG ID:lane1 PL:illumina LB:lib1 SM:sample1 genome.fasta clean_1.fq.gz clean_2.fq.gz | samtools view -S -b - sample1.bamID測序lane編號從fq文件名獲取如lane1PL測序平臺必須是GATK認可的如illumina、COMPLETE不能寫“CG”“MGI”SM樣本ID唯一標識多樣本分析時必用LB文庫名可選從測序報告獲取避坑重點-R參數(shù)的4個核心信息 平臺寫錯會報“not a recognized platform”錯誤后期改起來很麻煩命令解析-t 4用4個線程加速管道符|直接將比對結(jié)果SAM格式轉(zhuǎn)給Samtools用-b轉(zhuǎn)為二進制BAM格式節(jié)省空間后續(xù)分析更高效。2.2 數(shù)據(jù)排序按染色體位置“排好隊”BWA比對后的BAM文件是按read的測序順序排列的而后續(xù)分析需要按染色體位置排序。用Samtools完成命令samtools sort - 4 -m 4G -O bam -o sample1.sorted.bam sample1.bam參數(shù)說明-m 4G限制每個線程用4G內(nèi)存避免服務(wù)器內(nèi)存溢出文件名加“sorted”標識后續(xù)好區(qū)分。排序后文件會略小是壓縮算法導(dǎo)致的內(nèi)容無損失。2.3 標記重復(fù)序列剔除PCR擴增的“贗品”建庫時的PCR擴增會產(chǎn)生大量重復(fù)序列這些序列會干擾變異檢測增大假陽/假陰率必須標記或去除。主流用Picard的MarkDuplicates默認只標記不刪除更靈活。picard MarkDuplicates Isample1.sorted.bam Osample1.sorted.markdup.bam Msample1.markdup_metrics.txt參數(shù)說明I是輸入文件O是輸出文件M是重復(fù)序列統(tǒng)計報告可查看重復(fù)率一般低于30%算正常。如果非要刪除重復(fù)序列加REMOVE_DUPLICATEStrue參數(shù)即可。2.4 構(gòu)建索引讓工具“隨機訪問”文件標記重復(fù)后的BAM文件需要建索引方便后續(xù)工具快速定位特定區(qū)域。同時要給參考基因組做GATK專用索引兩步命令# 給BAM文件建索引 samtools index sample1.sorted.markdup.bam # 給參考基因組建GATK索引生成.dict和.fai文件 gatk CreateSequenceDictionary -R genome.fasta -O genome.dict samtools faidx genome.fasta運行后會生成sample1.sorted.markdup.bam.baiBAM索引、genome.dict和genome.fasta.fai參考基因組索引這三個文件缺一不可。2.5 變異檢測從GVCF到最終VCFGATK的HaplotypeCaller是目前最優(yōu)的變異檢測工具支持單樣本和多樣本分析核心分“生成GVCF→合并→基因型推斷”三步。步驟1單樣本生成GVCFGVCF文件包含所有位點信息無論是否變異便于后續(xù)多樣本合并分析。命令gatk HaplotypeCaller -R genome.fasta -I sample1.sorted.markdup.bam --emit-ref-confidence GVCF --min-base-quality-score 10 -O sample1.chr1.g.vcf.gz如果樣本多、染色體多建議寫shell腳本批量運行循環(huán)修改染色體號和樣本名效率翻倍。步驟2合并多樣本GVCF多個樣本按染色體合并先把同染色體的GVCF文件名存成列表再用CombineGVCFs合并# 生成GVCF列表 ls *.chr1.g.vcf.gz chr1_gvcf.list # 合并 gatk CombineGVCFs -R genome.fasta -V chr1_gvcf.list -L 1 -O chr1.merged.g.vcf.gz步驟3基因型推斷生成VCF將合并后的GVCF轉(zhuǎn)為最終的變異文件VCF包含SNP和InDel信息gatk GenotypeGVCFs -R genome.fasta -V chr1.merged.g.vcf.gz -O chr1.genotype.vcf.gz2.6 變異過濾剔除假陽性保留可靠結(jié)果剛生成的VCF是“原始數(shù)據(jù)”包含大量假陽性變異需要過濾。分SNP和InDel兩類處理非人類物種建議用“硬過濾”人類可用VQSR依賴已知變異集。# 提取SNP gatk SelectVariants -R genome.fasta -V chr1.genotype.vcf.gz -O chr1.snp.vcf -select-type SNP # 過濾SNP核心參數(shù)可根據(jù)數(shù)據(jù)調(diào)整 gatk VariantFiltration -V chr1.snp.vcf -O chr1.snp.filter.vcf -R genome.fasta --filter-expression QD 2.0 || FS 60.0 || MQ 40.0 --filter-name SNP_filter過濾參數(shù)說明QD變異質(zhì)量值、FS堿基偏倚、MQ比對質(zhì)量這些是GATK推薦的核心指標過濾后標記為“SNP_filter”的位點就是需要剔除的假陽性。3 收尾結(jié)果文件整理與后續(xù)分析方向上游分析結(jié)束后核心產(chǎn)出是過濾后的VCF文件如chr1.snp.filter.vcf后續(xù)可根據(jù)研究目的開展分析群體遺傳分析用PLINK做PCA、親緣關(guān)系分析用Admixture做群體結(jié)構(gòu)分析候選基因篩選結(jié)合注釋文件如ANNOVAR篩選位于外顯子區(qū)的有害變異關(guān)聯(lián)分析與表型數(shù)據(jù)結(jié)合做GWAS全基因組關(guān)聯(lián)分析定位性狀相關(guān)位點必看避坑總結(jié)軟件版本要匹配GATK4不兼容GATK3的命令安裝時明確指定版本conda install gatk44.4.0Read Group別瞎寫PL參數(shù)必須是GATK認可的SM參數(shù)要唯一否則后續(xù)報錯文件命名有規(guī)律建議用“樣本名_處理步驟.bam”格式如sample1_sorted.markdup.bam避免后續(xù)混淆服務(wù)器資源要算夠排序和變異檢測很耗內(nèi)存100G數(shù)據(jù)建議至少用16線程32G內(nèi)存中間文件別亂刪索引文件.bai/.fai和統(tǒng)計報告.metrics.txt后續(xù)可能用得上定期備份再清理