---
# System prepended metadata

title: WES analysis
tags: [VGHTC]

---

# 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