# Project 2: Prostate cancer WES analysis (Family-based study) ###### tags: `VGHTC` 從網頁版定序資料載 data ``` wget --user=THHsiaoLab --password=bEtAye7p --cut-dirs=2 -m -nH -np http://ngs.genomics.com.tw/download2/VGHT_THHsiaoLab/1709AHF-0001/ ``` ## Step 1: 事前準備 1. 將含有相同 sample name 的 R1/R2.fastq.gz 檔案依照 sample name 創建資料夾（Sample_xxx, PS: 會有 script 可以大量跑） 2. 將于晴的 script 全部複製一份到還有分析完 data 存到分析完的資料夾中（可以與 .fastq.gz 資料夾的地方不同） 3. 修改 1-1 的內的 input 路徑(有 Sample 資料夾的地方) 4. 修改 1-1 R1FILEPATTERN="R1fastq.gz", R2FILEPATTERN="R2fastq.gz", R1RMPATFORSAMPLEID="R1fastq.gz", R2RMPATFORSAMPLEID="R2fastq.gz" 與目標 fastq.gz 的命名方式相符 ## Step 2: run script 1-1 內包含： 1. fastQC 看 raw reads quality 2. Trim reads by trimmomatic 3. Alignment by BWA mem 4. Build index 5. run fastQC after QC steps 建議使用 screen 背景跑才跳出做其他事 ``` screen ``` 跑 1-1 (同時會跑 fastq_to_gvcf_part1_NBV_ENTV.sh) ``` sh 1-1.run_fastq_to_gvcf_part1_v3_NBV_ENTV_0710.sh ``` 看現在有幾個 screen ``` screen -ls screen -r xxxxx screen -D xxxxx #強制斷掉正在使用 screen 的使用者並重新連接 ``` 跳離 screen ``` ctrl+A+D ``` 關閉 screen ``` exit ``` gVCF 下各個 sample 的獨立資料夾中還會產生四個資料夾： 1. bwa 2. cleanfq 3. fastqc (200714增加的) 4. gatk 外層還會有會有 .log 檔會記錄所有跑完的資料（可以從中擷取 raw reads cleanup 的資訊）。 ### 詳細資料 1. fastqc 的資料夾內會有.html 的檔案可以用 FileZilla 傳到本機端用 FireFox 開啟。 - Basic Statistics 包含文件名(Filename)、文件類型、總序列數(Total Sequences)、序列長度(Sequence length)這些基本信息。 - Per base sequence quality 序列每個鹼基的平均質量，越高越好，需要注意會有部分序列在開頭部分質量差，所以根據這個圖在做質控時選擇兩端都去低質量或只去3’末端。 - Per tile sequence quality 新版本增加的功能，不太清楚是乾啥的。 - Per sequence quality scores 序列平均得分的數量。右側越高越好，也就是大多數序列質量都得分在30以上。均值低於27(也就是錯誤率大於0.002)記為警告，均值低於20(錯誤率0.01)記為不合格。 - Per base sequence content 顯示鹼基比例。正常情況下鹼基比例應該差不多。AT與CG差異大於10%記為警告，大於20%記為不合格。 - Per sequence GC content 每條序列GC含量百分比與模式化的正態分布GC含量相比較，超過15%記為警告，超過30%為不合格。 - Per base N content 每個鹼基位置N（未測到或不確定鹼基）的比例 - Sequence Length Distribution 各序列長度佔比 - Sequence Duplication Levels 序列重復級別 理想情況下所有序列應該是被隨機打斷了，測序後理論上不太會出現完全相同的序列。但由於PCR擴增或者污染等原因會造成一些重復序列，這裡顯示重復序列的比例。為了節省內存和時間，fastqc僅抽取了前100,000條序列，並只要超過75bp的序列就會被截斷到50bp來分析。藍線代表總reads的重復情況，紅色線代表Deduplicated reads的重復情況 一個好的結果應該是紅線藍線最左側越大越好。通常會在紅色線中看到一個高重復的峰，但同時在藍色線上消失，這表明重復序列並不顯著。如果在藍色的線中有峰值，說明在大量不同的高度重復序列，這可能是污染或嚴重的技術重復。如果非唯一序列數超過總數的50%就會被軟件判斷為不合格。由於測序技術限制，如果在這一步出現重復不及格，請注意在後續markduplicate的步驟去除重復序列。 - Overrepresented sequence 過表達序列。一般認為打斷的序列只有很少部分會重復，如果一段序列重復達到總值的0.1%記為警告，超過1%記為不合格。 2. trimmomatic 工具是用於illumina二代測序數據的reads處理，主要對接頭(adapter)序列和低質量序列進行過濾。 `https://www.twblogs.net/a/5d1129a0bd9eee1e5c820a54` - 參數說明 PE 設置使用trimmomatic處理雙端數據，單端數據用(‘SE’) -thread 16 設置線程數爲16 -phred33 設置鹼基的質量格式(默認-phred64，自v0.32版本之後可自動識別是phred33還是phred64) -trimlog trim.log 設置trimmommatic工具處理的日誌文件爲’trim.log’，每兩行爲一對reads信息 生成的日誌文件’trim.log’包含5列信息： 1) read名稱 2) 切除後剩餘的read長度 3) 切除後第一個鹼基所在位置(0-base)，也就是起始位置切除了幾個鹼基 4) 切除後最後一個鹼基所在位置 5) read末尾切除了幾個鹼基 PS由於該生成文件較大，如後續步驟不需該文件信息，可考慮不設置 * input_forward_R1.fq.gz 輸入的forward正向鏈R1對應的fastq文件 * input_reverse_R2.fq.gz 輸入的reverse反向鏈R2對應的fastq文件 * output_forward_paired_R1.fq.gz 處理後輸出的R1對應reads成對fastq文件 * output_forward_unpaired_R1.fq.gz 處理後輸出的R1對應reads不成對的fastq文件 * output_reverse_paired_R2.fq.gz 處理後輸出的R2對應reads成對fastq文件 * output_reverse_unpaired_R2.fq.gz 處理後輸出的R2對應reads不成對的fastq文件 * ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true 切除illumina接頭參數設置。說明ILLUMINACLIP的各參數之前，先解釋一下要使用的兩種切除接頭的模式和比對分值計算方法 * 下面是ILLUMINACLIP的各參數說明： ":<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>" ":TruSeq3-PE.fa" # <fastaWithAdaptersEtc>是fasta格式的接頭序列文件路徑 ":2" # <seed mismatches>是將接頭序列的一段(不超過16bp)作爲seed，與reads進行比對，能夠容忍的最大錯配(mismatch)數，這裏是最多2個錯配 ":30" # <palindrome clip threshold>是 a-R1, a-R2的比對分值閾值，達到閾值，才進行切除，這裏設置閾值爲30（大約比對50bp鹼基） ":10" # <simple clip threshold>是任意(接頭)序列與read比對最低分閾值，大於這個閾值，才進行切除，這裏設置閾值爲10（大約比對17bp鹼基） ":8" # <minAdapterLength>只作用於Palindrome模式，是設置檢測到接頭序列的最小長度（默認爲8，甚至可設置爲1） ":true" # <keepBothReads>只作用於palindrome模式，是設置是否保留反向鏈，這裏是說去除接頭序列後，由於正反鏈包含相同的序列信息(儘管序列是反向互補的)，默認情況(":false")下會去除反向鏈，設置爲":True"則保留反向鏈 * 參數option： -5, --seqType INT Sequence fq name type, 0->old fastq name, 1->new fastq name [默認0] Old fastq name:@FCD1PB1ACXX:4:1101:1799:2201#GAAGCACG/2 New fastq name: @HISEQ:310:C5MH9ANXX:1:1101:3517:2043 2:N:0:TCGGTCAC -1 SE的fastq文件 -2 PE的_2.fastq⽂件 （SE就不用輸-2） -T 核數，不同核數會有不同輸出結果，同⼀一套數據請設定同⼀一個核數 -l low quality threshold [5] -q low quality rate [0.5] -n N rate threshold [0.05] -f 5‘引物 -r 3’引物 （請填入已知的引物序列） -Q 2 設定sanger測序模式 -G 2 設定sanger測序模式，採用phred33 -o 輸出地址，按樣本來分⽂文件夾 -C 必填，輸出的fq1⽂件名(.gz格式) -D 輸出的fq2⽂件名(.gz格式) 3. Samtools是處理Alignment軟體產生出SAM檔之工具組 Usage: samtools <command> [options] Command: view SAM<->BAM conversion sort sort alignment file mpileup multi-way pileup depth compute the depth faidx index/extract FASTA tview text alignment viewer index index alignment idxstats BAM index stats (r595 or later) fixmate fix mate information flagstat simple stats calmd recalculate MD/NM tags and '=' bases merge merge sorted alignments rmdup remove PCR duplicates reheader replace BAM header cat concatenate BAMs targetcut cut fosmid regions (for fosmid pool only) phase phase heterozygotes ## Step 3: 準備 docker 環境 確定 docker 可以運作 ``` docker ``` 第一次使用 GATK 要先 pull image 故要修改 1-2 的註解位置，只留 docker pull broadinstitute/gatk:4.1.2.0 ``` sh 1-2.start_docker_for_gatk_NBV_ENTV.sh ``` 無權限時改用 "sudo" ``` sudo ./1-2.start_docker_for_gatk_NBV_ENTV.sh ``` 載完後再將 1-2 註解位置改回 docker run -v /data/yuching/Exome_data_20191219:/gatk/Exome_data_20191219 -it broadinstitute/gatk:4.1.2.0 並將路徑改為欲放入 container 的資料夾路徑 PS: 冒號後可改為方便自己閱讀的文字即可 ``` sudo ./1-2.start_docker_for_gatk_NBV_ENTV.sh ``` 進入 docker 之後把指令版本改成 bash，進入包含有 script 的資料夾，並 run 1-3，記得按 **NO** ``` bash sh 1-3.docker_switch_from_dash_to_bash.sh ``` ## Step 4: run GATK 1. 修改 1-4.run_fastq_to_gvcf_part2_NBV_ENTV.sh 中的 path 以及 R1 & R2 的命名方式 2. 確定 GATK 的 ref 在 docker 包含的資料夾下 3. 修改 fastq_to_gvcf_part1_NBV_ENTV.sh 中 ref 的 pathway 4. run 1-4.run_fastq_to_gvcf_part2_NBV_ENTV.sh ``` screen sh 1-4.run_fastq_to_gvcf_part2_NBV_ENTV.sh ``` # Additional process 從 Ensembl 下載 gtf 欓 ref: http://www.360doc.com/content/20/0827/13/19913717_932471349.shtml ref: http://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/ ``` http://ftp.ensembl.org/pub/grch37/release-101/gtf/homo_sapiens/ wget http://ftp.ensembl.org/pub/grch37/release-101/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.gtf.gz gunzip Homo_sapiens.GRCh37.87.gtf.gz ``` 從 ref.gtf 檔案中擷取部分並轉成 bed 欓 ``` vim my_gene_and_transcript_list.txt grep -f my_gene_and_transcript_list.txt Homo_sapiens.GRCh38.100.gtf > selected_genes.gtf awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' selected_genes.gtf | gtf2bed - > selected_genes.bed ``` 執行 runvcftools.sh 依照製作的依照製作的 bed 擷取我們資料內的位點 ``` sh runvcftools.sh ``` 若vcf非為合併欓，先列出 vcf 的檔案有哪些 ``` find . -iname "*.vcf" ``` 用 bgzip 將 vcf 轉為 vcf.gz 並列出來 ``` ls *.vcf | xargs -n1 bgzip find . -iname "*.vcf.gz" for F in *.vcf.gz ; do tabix -f -p vcf ${F} ; done bcftools merge Home/data/*vcf.gz -Oz -o Merged.vcf.gz ```