基因體 / 台大基因定序案 (NTUH Projects) === ###### tags: `基因體`, `生物資訊, `台大基因定序案` - NTUH: National Taiwan University Hospital <br> [TOC] <br> ## README ### 執行環境建置 ``` $ apt-get update $ apt-get install openjdk-8-jdk # optional $ update-alternatives --set java /usr/lib/jvm/java-8-openjdk-amd64/jre/bin/java $ apt-get install -y numactl graphviz parallel ``` <br> ### 認識測試資料 - SureSelection script + D15780 (測試的小資料) - SureSelection script + WES-EDA-013A (中資料) - WGS script + WGS-LIS-AI018A (大資料) ### 執行環境配置 1. 重新定義資料夾路徑，解壓縮放到 ```misc``` 下 ```everythings-misc.zip -> ./Everything/misc``` 2. 在 ```Everything``` 資料夾下，建立任意專案名稱，如 ```tj_project``` ```./Everything/tj_project``` 3. 在 ```tj_project``` 資料夾下，建立 reads，來放置測試資料 ``` └── Everything └── tj_project ├── my.ntuh.pipeline.b37.SureSelect.benchmark.haplotypecaller.spark.sh └── reads ├── D15780_S13_L001_R1.fastq.gz └── D15780_S13_L001_R2.fastq.gz ``` <br> 其中 xxx.gz 取自如下，選擇想要測試的資料； (可以先選擇小資料進行環境測試) ``` └── work ├── D15780_S13_L001_R1.fastq.gz (測試的小資料) <--- ├── D15780_S13_L001_R2.fastq.gz (測試的小資料) <--- ├── WES-EDA-013A_R1.fastq.gz (中資料) ├── WES-EDA-013A_R2.fastq.gz (中資料) ├── WGS-LIS-AI018A_R1.fastq.gz (大資料) └── WGS-LIS-AI018A_R2.fastq.gz (大資料) ``` <br> benchmark script 取自如下： ``` └── script ├── ntuh.pipeline.b37.SureSelect.benchmark.deepvariant.sh ├── ntuh.pipeline.b37.SureSelect.benchmark.haplotypecaller.spark.sh <--- ├── ntuh.pipeline.b37.WGS.benchmark.deepvariant.sh └── ntuh.pipeline.b37.WGS.benchmark.haplotypecaller.spark.sh ``` <br> 4. 回到 ```Everythings``` 目錄下，使用 soft link 建立絕對路徑，讓 script 可以跑 > 因為所寫的 script，是用絕對路徑寫死： > /Everythings/misc/nextflow-v18.10.1/nextflow run ... ``` Everythings$ sudo ln -s `pwd` /Everythings ``` <br> 5. 在剛剛建立的 ```tj_project``` 目錄下，執行 benchmark ``` Everythings/tj_project$ sh ntuh.pipeline.b37.SureSelect.benchmark.haplotypecaller.spark.sh & ``` bwa 會用掉 5.3GB java x 4 會用掉 1 GB <br> 6. 執行完，執行目錄下會有 report.html - 若成功，則沒有 error - terminal log 會顯示： ``` Success : true Exit status : 0 Error report: - ``` - report.html 會顯示： Nextflow workflow report Workflow execution completed successfully! - 各程序所花的時間 | process | duration | duration(%) | | -------- | -------- | -------- | | runBWA | 51.3s | 6.22% | | runSortSam | 1m 57s | 14.19% | | runMarkDuplicates | 1m 57s | 14.19% | | runBaseRecalibrator | 4m 49s | 35.06% | | runApplyBQSR | 1m 57s | 14.19% | | runHaplotypeCaller | 2m 13s | 16.13% | | TOTAL | 13m 44s | 100% | <br> ## Nextflow workflow report - Nextflow command ```bash nextflow run /Everythings/misc/config/ntuh.benchmark.haplotypecaller.spark.nf -c /Everythings/misc/config/config.nf --reads reads/*R{1,2}.fastq.gz --reference b37 --panel SURESELECT -w /tmp -with-dag flowchart.png ``` - each process - ### runBWA ``` numactl --interleave=all /Everythings/misc/bwa/bwa mem -M -R '@RG\tID:D15780_S13_L001\tSM:D15780_S13_L001\tPL:Illumina' -t 2 /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta D15780_S13_L001_R1.fastq.gz D15780_S13_L001_R2.fastq.gz | /Everythings/misc/samtools/samtools view -@ 2 -1 -o D15780_S13_L001.bam ``` <br> - ### runSortSam ``` /Everythings/misc/gatk-4.1.0.0/gatk SortSamSpark --input D15780_S13_L001.bam --output D15780_S13_L001.sorted.bam --sort-order coordinate --java-options "-XX:+UseNUMA -Xmx16G" --tmp-dir . -- --spark-runner LOCAL --spark-master local[4] --conf spark.local.dir=./tmp ``` - 指令說明 - [GATK4 command-line syntax](https://gatk.broadinstitute.org/hc/en-us/articles/360035531892-GATK4-command-line-syntax) <br> - ### runMarkDuplicates ([summary](https://gatk.broadinstitute.org/hc/en-us/articles/360036897012-MarkDuplicatesSpark-BETA-)) ``` /Everythings/misc/gatk-4.1.0.0/gatk MarkDuplicatesSpark --input D15780_S13_L001.sorted.bam --output D15780_S13_L001.sorted.dedup.bam --metrics-file D15780_S13_L001.markDuplicates.metrics --read-index D15780_S13_L001.sorted.bam.sbi --java-options "-XX:+UseNUMA" --tmp-dir . -- --spark-runner LOCAL --spark-master local[4] --conf spark.local.dir=./tmp ``` - 討論 - [Question: Why Do We Need Markduplicates For Variants Detection In Gatk Processing Pipeline? ](https://www.biostars.org/p/18784/) <br> - ### runBaseRecalibrator ([summary](https://gatk.broadinstitute.org/hc/en-us/articles/360036898312-BaseRecalibrator)) ``` /Everythings/misc/gatk-4.1.0.0/gatk BaseRecalibratorSpark --input D15780_S13_L001.sorted.dedup.bam --known-sites /Everythings/misc/bundle/b37/dbsnp_138.b37.vcf --known-sites /Everythings/misc/bundle/b37/1000G_phase1.indels.b37.vcf --known-sites /Everythings/misc/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf --output D15780_S13_L001.recal_data.grp --reference /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta --intervals /Everythings/misc/bed/b37/Agilent_SureSelect_Human_All_Exon_V7_S31285117_Regions.b37.bed --read-index D15780_S13_L001.sorted.dedup.bam.sbi --java-options "-XX:+UseNUMA" --tmp-dir . -- --spark-runner LOCAL --spark-master local[4] --conf spark.local.dir=./tmp ``` <br> argument details: ``` /Everythings/misc/gatk-4.1.0.0/gatk BaseRecalibratorSpark --input D15780_S13_L001.sorted.dedup.bam --known-sites /Everythings/misc/bundle/b37/dbsnp_138.b37.vcf --known-sites /Everythings/misc/bundle/b37/1000G_phase1.indels.b37.vcf --known-sites /Everythings/misc/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf --output D15780_S13_L001.recal_data.grp --reference /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta --intervals /Everythings/misc/bed/b37/Agilent_SureSelect_Human_All_Exon_V7_S31285117_Regions.b37.bed --read-index D15780_S13_L001.sorted.dedup.bam.sbi --java-options "-XX:+UseNUMA" --tmp-dir . -- --spark-runner LOCAL --spark-master local[4] --conf spark.local.dir=./tmp ``` - ```--intervals```: - [Interval-related arguments and syntax](https://software.broadinstitute.org/gatk/documentation/article?id=11009) - ```--java-options```: - [4. Adding Java arguments](https://software.broadinstitute.org/gatk/documentation/article?id=11050) - [Java8 JVM参数解读](https://www.zybuluo.com/changedi/note/975529) - empty - 15.0 min - ```-Xmn2g``` - 15.82 min - ```-Xmn2g -Xms2g``` - 13.68, 14.43, 15.05, 13.74, 14.52, 13.43 min - ```-XX:+UseNUMA``` >使用NUMA[5]开启性能优化。默认不开启，该项只有在开启了-XX:+UseParallelGC后才有效。 - 15.11, 15.42, 14.92 min - ```-XX:+AggressiveHeap -XX:+UseParallelGC -XX:+UseNUMA``` - 15.27 min - ```-XX:+AggressiveHeap -XX:+UseParallelGC -XX:+UseNUMA``` - 16.33, 13.97, 15.90 min <br> 單獨執行指令時，若遇到底下錯誤，請把檔案路徑改成絕對路徑 > A USER ERROR has occurred: Failed to read bam header from D15780_S13_L001.sorted.dedup.bam Caused by:File D15780_S13_L001.sorted.dedup.bam does not exist - **```--input```** ``` --input D15780_S13_L001.sorted.dedup.bam ``` 改成(舉例) ``` --input /home/tj/Asus/workplace/Gene/Everythings/myproject/OUTPUT/1.ALIGNMENT/D15780_S13_L001.sorted.dedup.bam ``` - **```--read-index```** ``` --read-index D15780_S13_L001.sorted.dedup.bam.sbi ``` 改成(舉例) ``` --read-index /home/tj/Asus/workplace/Gene/Everythings/myproject/OUTPUT/1.ALIGNMENT/D15780_S13_L001.sorted.dedup.bam.sbi ``` - **testing reports:** - Elapsed time: 14.12, 13.67, 14.07 (min) - CPU: low <br> - ### runApplyBQSR ([summary](https://gatk.broadinstitute.org/hc/en-us/articles/360037268511-ApplyBQSR)) BQSR: Base Quality Score Recalibration ``` /Everythings/misc/gatk-4.1.0.0/gatk ApplyBQSRSpark --bqsr-recal-file D15780_S13_L001.recal_data.grp --input D15780_S13_L001.sorted.dedup.bam --output D15780_S13_L001.sorted.dedup.recal.bam --intervals /Everythings/misc/bed/b37/Agilent_SureSelect_Human_All_Exon_V7_S31285117_Regions.b37.bed --reference /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta --read-index D15780_S13_L001.sorted.dedup.bam.sbi --java-options "-XX:+UseNUMA" --tmp-dir . -- --spark-runner LOCAL --spark-master local[4] --conf spark.local.dir=./tmp ``` <br> - ### runHaplotypeCaller ``` if grep -q "avx" /proc/cpuinfo; then HMM=AVX_LOGLESS_CACHING_OMP else HMM=LOGLESS_CACHING fi if grep -q "avx2" /proc/cpuinfo; then SW=AVX_ENABLED else SW=JAVA fi numactl --interleave=all /Everythings/misc/gatk-4.1.0.0/gatk HaplotypeCaller -R /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta -I D15780_S13_L001.sorted.dedup.recal.bam -L /Everythings/misc/bed/b37/Agilent_SureSelect_Human_All_Exon_V7_S31285117_Regions.b37.bed --dbsnp /Everythings/misc/bundle/b37/dbsnp_138.b37.vcf -O D15780_S13_L001.vcf --read-index D15780_S13_L001.sorted.dedup.recal.bam.bai --tmp-dir . --smith-waterman $SW --pair-hmm-implementation $HMM --native-pair-hmm-threads 2 ``` <br> argument details: ``` numactl --interleave=all /Everythings/misc/gatk-4.1.0.0/gatk HaplotypeCaller -R /Everythings/misc/bundle/b37/human_g1k_v37_decoy.fasta -I D15780_S13_L001.sorted.dedup.recal.bam -L /Everythings/misc/bed/b37/Agilent_SureSelect_Human_All_Exon_V7_S31285117_Regions.b37.bed --dbsnp /Everythings/misc/bundle/b37/dbsnp_138.b37.vcf -O D15780_S13_L001.vcf --read-index D15780_S13_L001.sorted.dedup.recal.bam.bai --tmp-dir . --smith-waterman $SW --pair-hmm-implementation $HMM --native-pair-hmm-threads 2 ``` <br> <br> ## 從 ESC4000(10.78.26.241) 下載資料的速度 ``` $ scp diatango_lin@10.78.26.241:/Everythings/dataset/* . D15780_S13_L001_R1.fastq.gz 100% 58MB 3.4MB/s 00:17 D15780_S13_L001_R2.fastq.gz 100% 64MB 3.4MB/s 00:19 WES-EDA-013A_R1.fastq.gz 100% 2536MB 3.2MB/s 13:06 WES-EDA-013A_R2.fastq.gz 100% 2751MB 3.2MB/s 14:08 WGS-LIS-AI018A_R1.fastq.gz 100% 33GB 3.3MB/s 2:50:09 WGS-LIS-AI018A_R2.fastq.gz 100% 34GB 3.3MB/s 2:56:59 ``` ``` ├── [ 60858292] D15780_S13_L001_R1.fastq.gz ├── [ 66890769] D15780_S13_L001_R2.fastq.gz ├── [ 2658993692] WES-EDA-013A_R1.fastq.gz ├── [ 2884234372] WES-EDA-013A_R2.fastq.gz ├── [35314302136] WGS-LIS-AI018A_R1.fastq.gz └── [36820690257] WGS-LIS-AI018A_R2.fastq.gz ``` - 檔案大小，全部共 77805969518 bytes (=72.46 GB) - 下載時間，全部共 6h:14m:58s - 內網平均下載速度： 77805969518 / ((6*60+14)*60+58) / 1024 / 1024 = 3.2981 (MB/s)