WES analysis - HackMD

# WES analysis ###### tags: `VGHTC` for THRV-TOPMed program bcf file PWD: ``` ssh @172.22.116.45 root@joinet cd /data3/Project/President_Hsu_Taichi/THRV/freeze.8/minDP0 ``` ## 0. upload fastq to basespace intro: https://hackmd.io/z3iV2Oa-RCmzNHAsGgomKw https://developer.basespace.illumina.com/docs/content/documentation/cli/cli-examples#Install Jack workbook: https://hackmd.io/uD9RiX5ZRp6SvhxXVv36Dw?view Jack workbook: https://hackmd.io/VLkjxKzCQJSZDJf_tPyfMA illumina workbook: https://developer.basespace.illumina.com/docs/content/documentation/cli/cli-examples#Listprojects ### 0-1. make new project in basespace 1. log in to basespact and creat new project 2. move to the server to print out the project list ``` bs list projects bs list datasets ``` 3. rename R1 R2 file name match to basespace Make file names (with path) in a text file. ``` pwd | xargs -I % find % -type f -iname '*.fastq.gz'> ../filename.txt ``` Copy to local and add a new column with new names. PS: drop the file type label in new name col (ex: old.txt > new) Upload the rename file. ``` awk '{f=$1;sub(".*\\.",$2".",$1);system("mv "f" "$1)}' new_filename3.txt ``` 4. list the files in a text file and copy to another folder Jupyter notebook: /home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/Find file path and size manually delete size small duplicates ``` while read file; do cp "$file" /data3/prostate_cancer/germline/rawData_for_BaseSpace; done < /data3/prostate_cancer/germline/cleaned_data.txt ``` ``` while read file; do mv "$file" /data3/prostate_cancer/germline/reload_rawData_for_BaseSpace; done < /data3/prostate_cancer/germline/cleaned_data_210928.txt ``` ``` while read file; do rm -Rf "$file" ; done < done_list.txt ``` ### 0-2 upload fastq.gz files ``` bs list projects bs upload dataset -p 287706419 --recursive . ``` ### 0-3 check sample list ``` bs list projects bs list datasets --project-id=287706419 ``` Delete wrong files ``` bs delete dataset -i ds.835b7f9f4bf14c42aac120f6a89555ca ``` if more then 10 do for loop in /data3/prostate_cancer/germline/del.sh ### 0-4 find biosample ID ``` bs biosample list --project-id=287706419 ``` ### 0-5 do analysis <WES DRAGEN Germline> https://support.illumina.com/content/dam/illumina-support/help/Illumina_DRAGEN_Bio_IT_Platform_v3_7_1000000141465/Content/SW/Informatics/Dragen/GPipelineVarCalFilt_fDG.htm 列出 DRAGEN 分析的所有可設定參數 ``` bs application launch --name "DRAGEN Germline" --app-version "3.8.4" --list ``` * target_bed_id: 提供一個 list 指定要分析的位點 * coverage_list, coverage_list.XXX:舊有的QC選項，目前已合併到 multiqc_checkbox 功能中 * cnv_XXX: copy number variation analysis (若要做就把 checkbox 參數改成 1) * sv_XXX: structure varuantion analysis (若要做就把 checkbox 參數改成 1) * sv-force: 強迫分析某一區段 * eh_checkbox: Expansion Hunter 針對 repeat 做計算 * dupmark_checkbox: Duplicate Marking 會 capture PCR 造成的 bias 並在後續分析去除 (若勾選，在 PCR-base 的 library 最後會沒東西；在 capture-base 的才能使用) * vc_XXX: 未知功能 * cyp2d6: 單獨針對 CYP2D6 gene 做 variant calling * metrics_checkbox: 各式計算結果 * nirvana: 由 illumina 設計的 variant annotation 程式 (https://illumina.github.io/NirvanaDocumentation/) * md5_all: 使用於資料傳輸前後的資料完整度確認 (會產生6-10個碼的 barcod) * dbsnp: given a list of SNPid to force variant calling * forcegt: given a list to force variant calling 無論有沒有 sample 有 mutation 都會列出來 * concordance_vcf: 可案照 DRAGEN 的格式更改 variant calling 的參數 * vc-filter: 選擇 vcf 檔案中 FILTER 的欄位是否篩選掉 preQC 不好的點 ``` bs application launch --name "DRAGEN Germline" \ --app-version "3.8.4" \ -o app-session-name:"210929-PC_anaset6" \ -o project-id:287706419 \ -o ht-ref:"hg38-altaware-cnv-anchor.v8" \ -o input_list.sample-id:437926492,437908471,437879444,437884447,437868431,437858422,437838401,437800363,437776347,437748321,437744307,437728292,437726289,437711275,437696260,437687250,437662225,437658221,437650215 \ -o sv_checkbox:1 \ -o eh_checkbox:1 \ -o dupmark_checkbox:1 \ -o cyp2d6:1 \ -o nirvana:1 ``` ### 0-6 download vcf files in a project ``` bs project download -o /data3/prostate_cancer/DRAGEN_results -n 210802_Prostate_cancer --extension=vcf.gz ``` list out the output files ``` pwd | xargs -I % find % -type f -iname '*.hard-filtered.vcf.gz'> ../211006_prostate_vcf.txt ``` extract less frequent and rare variants * run shell /home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/1.runbcftools_extract_position.sh * merge vcf.gz /home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/2.runbcftools_merge.sh Merge files ``` bcftools merge *vcf.gz -Oz -o ../DRAGEN_PC_minor0.05.vcf bcftools norm -m -both DRAGEN_PC_minor0.05.vcf > DRAGEN_PC_minor0.05.norm.vcf bgzip -c DRAGEN_PC_minor0.05.norm.vcf > DRAGEN_PC_minor0.05.norm.vcf.gz tabix -f DRAGEN_PC_minor0.05.norm.vcf.gz ``` Merge with control files ``` bcftools merge --m none noPC_monir0.05.vcf.gz DRAGEN_PC_minor0.05.norm.vcf.gz -Oz -o Minor0.05_merged.vcf.gz gunzip Minor0.05_merged.vcf.gz ``` check sex with plink ``` plink --vcf freeze.8.chrX.pass_and_fail.gtonly.minDP0.vcf.gz --allow-extra-chr --check-sex ``` ### 0-7. List out the samples ``` for sample in `bcftools view -h Minor0.05_merged_male_ACMG_filtered.recode.vcf | grep "^#CHROM" | cut -f10-`; do echo $sample; done > SampleNames.txt ``` ## 1. from bcf file to vcf file ``` bcftools view -Oz -o /home/amy/TPMI/Ori_Apt_Adv/WGS_count/bcf_file/freeze.8.chr19.pass_and_fail.gtonly.minDP0.vcf.gz /data3/Project/President_Hsu_Taichi/THRV/freeze.8/minDP0/freeze.8.chr19.pass_and_fail.gtonly.minDP0.bcf ``` bcf 依照 chromosome 拆分，所有的都已轉成 vcf file： ``` cd /data3/Project/President_Hsu_Taichi/THRV/freeze.8/vcf_gz ``` ## 2. Extract variants ### 2-1-1. prepare position file * No header * Two column * split with \t ``` chr1 123456 chr2 567890 chrX 147852 ``` ### 2-1-2. Prepare gene region BED file DOWNLOAD gene region (.gtf) from UCSC ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/ * /home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/Homo_sapiens.GRCh37.87.chr.gtf * /home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/Homo_sapiens.GRCh38.104.chr.gtf Run /home/amy/TPMI/Ori_Apt_Adv/Generate BED file.ipynb ``` #!/usr/bin/env python """ GTF.py Kamil Slowikowski December 24, 2013 Read GFF/GTF files. Works with gzip compressed files and pandas. http://useast.ensembl.org/info/website/upload/gff.html LICENSE This is free and unencumbered software released into the public domain. Anyone is free to copy, modify, publish, use, compile, sell, or distribute this software, either in source code form or as a compiled binary, for any purpose, commercial or non-commercial, and by any means. In jurisdictions that recognize copyright laws, the author or authors of this software dedicate any and all copyright interest in the software to the public domain. We make this dedication for the benefit of the public at large and to the detriment of our heirs and successors. We intend this dedication to be an overt act of relinquishment in perpetuity of all present and future rights to this software under copyright law. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. For more information, please refer to <http://unlicense.org/> """ from collections import defaultdict import gzip import pandas as pd import re GTF_HEADER = ['seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame'] R_SEMICOLON = re.compile(r'\s*;\s*') R_COMMA = re.compile(r'\s*,\s*') R_KEYVALUE = re.compile(r'(\s+|\s*=\s*)') def dataframe(filename): """Open an optionally gzipped GTF file and return a pandas.DataFrame. """ # Each column is a list stored as a value in this dict. result = defaultdict(list) for i, line in enumerate(lines(filename)): for key in line.keys(): # This key has not been seen yet, so set it to None for all # previous lines. if key not in result: result[key] = [None] * i # Ensure this row has some value for each column. for key in result.keys(): result[key].append(line.get(key, None)) return pd.DataFrame(result) def lines(filename): """Open an optionally gzipped GTF file and generate a dict for each line. """ fn_open = gzip.open if filename.endswith('.gz') else open with fn_open(filename) as fh: for line in fh: if line.startswith('#'): continue else: yield parse(line) def parse(line): """Parse a single GTF line and return a dict. """ result = {} fields = line.rstrip().split('\t') for i, col in enumerate(GTF_HEADER): result[col] = _get_value(fields[i]) # INFO field consists of "key1=value;key2=value;...". infos = [x for x in re.split(R_SEMICOLON, fields[8]) if x.strip()] for i, info in enumerate(infos, 1): # It should be key="value". try: key, _, value = re.split(R_KEYVALUE, info, 1) # But sometimes it is just "value". except ValueError: key = 'INFO{}'.format(i) value = info # Ignore the field if there is no value. if value: result[key] = _get_value(value) return result def _get_value(value): if not value: return None # Strip double and single quotes. value = value.strip('"\'') # Return a list if the value has a comma. if ',' in value: value = re.split(R_COMMA, value) # These values are equivalent to None. elif value in ['', '.', 'NA']: return None return value ``` ``` gene_file=dataframe('/home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/Homo_sapiens.GRCh37.87.chr.gtf') gene_all=pd.read_csv('/home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/Gene_list_cancer+DNArepair.txt',sep='\t') gene_list=gene_all.Gene.unique() gene_file_extract=gene_file[gene_file['gene_name'].isin(gene_list)] # print(gene_file_extract.shape) display(gene_file_extract.head()) gene_range=gene_file_extract[gene_file_extract['feature']=='gene'] gene_range['seqname']=gene_range['seqname'].astype(object) # gene_range['seqname']=gene_range.seqname.astype(str) gene_range['chrom'] = 'chr' + gene_range['seqname'].astype(str) # gene_range['chrom'] = gene_range['seqname'].apply(lambda x: 'chr'+str(x)) gene_range['chromStart'] = gene_range['start'] gene_range['chromEnd'] = gene_range['end'] gene_range['note'] = 'gene_id \"'+gene_range['gene_id'].astype(str)+'\"; exon_number \"'+gene_range['exon_number'].astype(str)+'\"; exon_id \"'+gene_range['exon_id'].astype(str)+'\";' #bed_file.to_csv('/home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/cancer_hallmarks.bed',header=False) #display(bed_file) gene_range[['chrom','chromStart','chromEnd','strand','gene_name']].to_csv('/home/amy/TPMI/Ori_Apt_Adv/Prostate_WES/DNAreapir_Cancer_gene_37.bed', header=False,index=False,sep='\t') ``` ### 2-2. Extract sites with position file and count variants ``` vcftools --gzvcf /data3/Project/President_Hsu_Taichi/THRV/freeze.8/vcf_gz/ \ --positions /home/amy/TPMI/Ori_Apt_Adv/WGS_count/position_all.txt \ --counts \ --out /home/amy/TPMI/Ori_Apt_Adv/WGS_count/recode_new/freeze.8.chr19.pass_and_fail.gtonly.minDP0 ``` filtered ``` bcftools view -i 'AN>=4340' Minor0.05_merged_ACMG.recode.vcf > Minor0.05_merged_ACMG_filtered.recode.vcf bcftools view -i 'AN>=2000' Minor0.05_merged_male_ACMG.recode.vcf > Minor0.05_merged_male_ACMG_filtered.recode.vcf ``` #### 2-2-1. Modify the annotation header of vcf file ``` grep -v "^##" freeze.8.chr19.pass_and_fail.gtonly.minDP0.recode.vcf > freeze.8.chr19.pass_and_fail.gtonly.minDP0.recode_1.vcf ``` #### 2-2-2. Generate vcf.gz file for merging ``` bgzip -c freeze.8.chr19.pass_and_fail.gtonly.minDP0.recode_1.vcf >freeze.8.chr19.pass_and_fail.gtonly.minDP0.recode_1.vcf.gz ``` #### 2-2-3. Generate vcf.gz.tbi for merging ``` tabix -f freeze.8.chr19.pass_and_fail.gtonly.minDP0.recode_1.vcf.gz ``` #### 2-2-4. merge vcf files https://shicheng-guo.github.io/bioinformatics/1923/02/28/bcftools-merge ``` vcf-merge /home/amy/TPMI/Ori_Apt_Adv/WGS_count/bcf_file/*recode_1.vcf.gz > /home/amy/TPMI/Ori_Apt_Adv/WGS_count/merge.vcf find ~+ -maxdepth 1 -type f -name '*.vcf' > ../cel_files.txt bcftools merge -l merge.txt -Oz -o merge.vcf.gz for sample in `bcftools view -h TWB_noPC_minor0.05.vcf | grep "^#CHROM" | cut -f10-`; do echo $sample; done > sampleNames_final.txt ``` #### 2-2-5. Calculate variant count with vcftools ``` vcftools --vcf merge.vcf --counts --out merge_counts ``` ## 3. Extract sites in a position range ``` vcftools --gzvcf /data3/Project/President_Hsu_Taichi/THRV/freeze.8/vcf_gz/freeze.8.chr19.pass_and_fail.gtonly.minDP0.vcf.gz \ --chr chr19 --from-bp 10950000 --to-bp 11100000 \ --recode-INFO-all --recode \ --out /home/amy/TPMI/Ori_Apt_Adv/WGS_count/Freq_analysis/LDLR_analysis ``` ## 4. Copy files from 172.22.116.45 to 172.22.116.52 ``` scp /home/amy/TPMI/Ori_Apt_Adv/WGS_count/Freq_analysis/STAP1_analysi.recode.vcf root@172.22.116.52:/data/amy/210316_FH/ rsync -avzh /home/amy/TPMI/Ori_Apt_Adv/WGS_count/Freq_analysis/STAP1_analysi.recode.vcf /root@172.22.116.52:/data/amy/210316_FH/ ``` ## 5. ANNOVAR ``` ssh root@172.22.116.52 1qaz@WSX pwd=/data/amy/ sh 1.0_make_avinput.sh sh 2.0_ANNOVAR.sh ``` * copy the output back to 45 * pwd=/home/amy/TPMI/Ori_Apt_Adv/WGS_count * run 1.Count_variants_in_WGS.ipynb ### 5-1. Convert .vcf file to .avinput file ``` perl /home/2019GenomicsEpidemiologyWorkshop/annovar/convert2annovar.pl -format vcf4 /data/amy/210316_FH/APOB_analysis.recode.vcf > /data/amy/210316_FH/ANNOVAR_input_file/APOB_analysis.recode.avinput ``` ### 5-2. Append info from databases ``` perl /data/amy/210316_FH/ANNOVAR_input_file/annovar_SCH.pl -t 8 -i /data/amy/210316_FH/ANNOVAR_input_file/APOB_analysis.recode.avinput -r hg38 -o /data/amy/210316_FH/ANNOVAR_results/recode_annovar ``` ## 6. 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 並列出來 merge https://shicheng-guo.github.io/bioinformatics/1923/02/28/bcftools-merge ``` ls *.vcf | xargs -n1 bgzip find . -iname "*.vcf.gz" for F in *.vcf.gz ; do tabix -f -p vcf ${F} ; done ``` 每個 sample 獨立 vcf ``` bcftools merge *vcf.gz -Oz -o ../Merged.vcf.gz ``` 多個 sample 分成不同的 vcf ``` bcftools concat *vcf.gz -Oz -o ../noPC_monir0.05.vcf.gz ``` 多個 sample 多種 vcf ``` bcftools merge --m none noPC_monir0.05.vcf.gz DRAGEN_PC_minor0.05.norm.vcf.gz -Oz -o Minor0.05_merged.vcf.gz ``` 活用 bcftools 1. normalized 將同一個 position 的不同 allele 紀錄分開 question https://www.biostars.org/p/210141/ workbook https://samtools.github.io/bcftools/howtos/query.html workbook https://samtools.github.io/bcftools/bcftools.html#norm ``` bcftools norm -m -both PC_minor0.05.vcf > PC_minor0.05_nor2.vcf ``` 2.Batch effect issues after merge of two datasets https://www.biostars.org/p/496434/ 3. 重新計算 Minor allele according to 我們自己的 population ### 6-1. SNPnexus annotation (生成 vcf file) https://www.snp-nexus.org/v4/ https://gwaslab.com/2021/07/14/rsid-chrpos-conversion/ ## 7. RENOVO https://github.com/mazzalab-ieo/renovo PS: 要跟 ANNOVAR 裝在同一個 server (172.22.116.52) ### 7-1 安裝 ``` git clone https://github.com/mazzalab-ieo/renovo.git cd renovo conda env create -f ReNOVo.yml conda activate renovo conda install -c r r-curl r-httr r-rvest r-readxl r-tidyverse conda install -c bioconda r-openxlsx python -m pip install scikit-learn==0.20.3 pip install pandas pip install matplotlib pip install seaborn pip install argparse Rscript -e "install.packages(c('openxlsx','tidyverse','readxl'), repos='http://cran.us.r-project.org')" sudo yum install -y libxml2-dev libcurl4-openssl-dev libssl-dev Rscript -e "install.packages(c('rvest','xml2'), repos='http://cran.us.r-project.org')" ``` ### 7-2 Change the interpreter (python/Rscript) PS: which APP (find the absolute path) in ReNOVo.py: #!/usr/bin/env python in preprocessing.R: #!/usr/bin/Rscript in Renovo_implementation.py:#!/usr/bin/env python change the script pathway * in Renovo_implementation.py rf = joblib.load("/data/amy/renovo/Files/RF_model.pkl") keep = pd.read_csv("/data/amy/renovo/Files/variables.txt",sep="\t") col_fin = pd.read_csv("/data/amy/renovo/Files/ordered_cols.txt",sep="\t") * in preprocessing.R source("/data/amy/renovo/Scripts/FixData_median.R") median<- read_excel("/data/amy/renovo/Files/median_correct.xlsx") PS: trouble in loading tidyverse ``` locate libicui18n.so.58 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/anaconda2/envs/renovo/lib/ ``` ### 7-3 Run RENOVO ``` conda activate renovo python /data/amy/renovo/ReNOVo.py -a /home/2019GenomicsEpidemiologyWorkshop -p /data/amy/211004_Prostate_cancer/ANNOVAR_input_file/ ``` ## 8 Variant Effect Predictor (VEP) ### 8-1 install https://github.com/Ensembl/ensembl-vep http://asia.ensembl.org/info/docs/tools/vep/script/vep_download.html ``` conda install -c bioconda ensembl-vep ``` ``` perl -version. cpan Archive::Zip cpan DBI conda install -c bioconda mysqlclient conda install -c conda-forge mysql which mysql export PATH=$PATH:/opt/conda/bin/mysql cpan DBD::mysql cd ensembl-vep Module::Build module conda install -c bioconda perl-bio-db-hts cpan Bio::DB::HTS ``` http://www.htslib.org/download/ https://www.biostars.org/p/328831/ ``` wget https://github.com/samtools/htslib/releases/download/1.9/htslib-1.9.tar.bz2 \ && tar -vxjf htslib-1.9.tar.bz2 \ && cd htslib-1.9 \ && make ``` PKG_CONFIG_PATH=$PKG_CONFIG_PATH:/usr/bin/pkg-config export PKG_CONFIG_PATH

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.