KoGES imputation (2) - QC 및 filtering #1

 

Imputation은 아래 내용 참조하여 시행하셨을 겁니다.

 

https://nephro.tistory.com/11

KoGES imputation (1)

 

 

이제 imputed 라는 디렉토리 안에 imputation 된 파일들이 존재할 것입니다. 예를 들어 log 파일 하나를 보도록 하겠습니다. 저는 root ID로 작업을 하도록 bash file을 만들었는데, nas에서 작업이 이루어지기 때문입니다. sudo su로 root로 변경한 후 작업을 해야, nas에서도 권한 문제로 귀찮아지지 않는 것 같습니다.

root@computerID:/home/samelsims/nas/Gene_data/KoGES/imputation/KCHIP_ver3/imputed# less Kchip_RSID_chr1_1K_imputed_beagle_b37_EAS.log

root@computerID (computerID는 임의로 하였습니다) 뒤에 현재 디렉토리가 나오고 #로 끝이 납니다. #로 끝난다는 것이 곧 root (즉, sudo)를 의미합니다.

 

 

명령어는 아래와 같이 log 파일을 보는 것입니다. 앞으로는 이렇게 명령어만 입력하도록 하겠습니다.

less Kchip_RSID_chr1_1K_imputed_beagle_b37_EAS.log

 

 

어쨌든 log 파일을 보면 아래와 같습니다.

 

 

beagle.05May22.33a.jar (version 5.4)
Copyright (C) 2014-2022 Brian L. Browning
Enter "java -jar beagle.05May22.33a.jar" to list command line argument
Start time: 01:04 PM KST on 16 Nov 2023

Command line: java -Xmx512000m -jar beagle.05May22.33a.jar
  gt=prephased/Kchip_RSID_chr1_1K_beagle_b37_EAS.vcf.gz
  ref=/home/samelsims/driveD/ref_gene/1K_phase3_b37_vcf/EAS/chr1.1kg.phase3.v5a.EAS.vcf.gz
  map=/home/samelsims/driveD/ref_gene/1K_phase3_b37_bref3/plink.chr1.GRCh37.map
  chrom=1
  out=imputed/Kchip_RSID_chr1_1K_imputed_beagle_b37_EAS
  impute=true
  window=40.0
  nthreads=18

Reference samples:                  504
Study     samples:               72,292

Window 1 [1:10177-18871042]
Reference markers:              225,316
Study     markers:                3,891

Imputation time:               31 minutes 33 seconds

Window 2 [1:18268831-53836636]
Reference markers:              362,442
Study     markers:                5,222

Imputation time:               33 minutes 50 seconds

Window 3 [1:49070164-87760656]
Reference markers:              406,182
Study     markers:                5,319

Imputation time:               45 minutes 18 seconds

Window 4 [1:85578740-144962195]
Reference markers:              384,313
Study     markers:                4,766

Imputation time:               2 hours 53 minutes 3 seconds

Window 5 [1:144224289-175818910]
Reference markers:              307,578
Study     markers:                4,243

Imputation time:               1 hour 46 minutes 3 seconds

Window 6 [1:173092350-213583338]
Reference markers:              431,578
Study     markers:                5,298

Imputation time:               1 hour 9 minutes 26 seconds

Window 7 [1:211080189-240239399]
Reference markers:              327,882
Study     markers:                5,146

Imputation time:               52 minutes 52 seconds

Window 8 [1:239251841-249240543]
Reference markers:              124,957
Study     markers:                1,798

Imputation time:               18 minutes 58 seconds

Cumulative Statistics:

Reference markers:            2,570,248
Study     markers:               35,683

Imputation time:               8 hours 51 minutes 3 seconds
Total time:                    8 hours 55 minutes 17 seconds

End time: 09:59 PM KST on 16 Nov 2023
beagle.05May22.33a.jar finished

 

 

보시면 chromosome 1 imputation 하는 데에만 9시간 가까이 걸린 것을 알 수 있습니다. nthreads 옵션은 사용한 쓰레드 수를 나타내는데, 제 컴퓨터가 I9, process 10개 짜리라 thread는 20개입니다. 18개 thread로 imputation을 했는데도 chromosome 1개 imputation하는데 그 정도 시간이 걸렸습니다. 

 

이렇게 chromosome 1 부터 22 까지 파일들이 만들어 지게 됩니다. 이후에 imputation 된 genotype data중에서 필요한 것들만 추려내야 하는데, 그 과정은 다음과 같습니다

 

(1) Tabix - 앞으로 있을 작업을 위한 색인 과정

(2) dosage  filtering - imputation accuracy (DR)과 allele frequency가 cut-off 이상의 정보만 추려내기

(3) duplication 제거 - duplication check와 check된 duplication을 제거하는 두 단계로 나뉩니다.

(4) merge - duplication까지 제거된 파일들을 하나로 합쳐주는 과정.

 

우선 오늘은 Tabix로 색인 과정과 dosage filtering 과정을 해보도록 하겠습니다.

 

1. Tabix

이후에는 indexing 작업이 필요합니다. 이는 향후 작업을 속도를 높이기 위해, 색인을 만드는 작업입니다. 이 작업을 위해서는 tabix라는 프로그램이 필요합니다. samtool이라는 툴을 만든 곳에서 같이 만들었습니다.

 

https://www.htslib.org/doc/tabix.html

tabix(1) manual page

 

우선 코드는 아래와 같습니다. 파일 명은 저의 경우에는 step5_1_tabix.sh 이런식으로 만들었습니다. imputation 과정 중에 여러가지 QC등을 하게 되는데, 그렇게 하다가 step4에서 imputation을 실제로 시행하고, step5_1 번째 QC 작업인 tabix를 실행한다는 의미를 주려구요.

#!/bin/bash

# name of file
FILE_NAME=Kchip_RSID

# output directory
mkdir imputed

# baseline setting
CHR=1
total_chr=22
thread=16

##########################################
# for in. thread control
for process in $(seq 1 $total_chr)
    do j_count=`jobs -l|wc|awk '{print $1}'`
    if [[ $j_count -ge $thread ]];then
        until [[ $j_count -lt $thread ]]
            do j_count=`jobs -l|wc|awk '{print $1}'`
            sleep 10
            done
    fi
##########################################

# Actual process of threads

/usr/bin/tabix imputed/${FILE_NAME}_chr${CHR}_1K_imputed_beagle_b37_EAS.vcf.gz &
echo "Chromosome ${CHR} is being processed using Tabix."

CHR=$(($CHR+1))

###########################################
    done

lastPIDs=`jobs -l|awk '{print $2}'`
wait $lastPIDs
echo "";echo "Indexing process have been done."

 

 

#!/bin/bash

# bash 임을 선언합니다. sh 파일 가장 맨 처음에 필요하죠.

# name of file
FILE_NAME=Kchip_RSID

# 뒤에서 사용할 파일명으로 여기서 미리 선언해 줍니다.

 

# output directory
mkdir imputed

# imputed라는 디렉토리를 만들어서, 이 디렉토리에 만들어진 파일들을 저장하도록 할 예정입니다.

## baseline setting
CHR=1
total_chr=22
thread=16
# 우선 시작 CHR을 1로 설정합니다. 그리고 총 chromosome 갯수를 22개로 설정합니다.
# thread는 16개로 설정하도록 하겠습니다.

##########################################
# for in. thread control
for process in $(seq 1 $total_chr)
  do j_count=`jobs -l|wc|awk '{print $1}'`
  if [[ $j_count -ge $thread ]];then
    until [[ $j_count -lt $thread ]]
      do j_count=`jobs -l|wc|awk '{print $1}'`
      sleep 10
      done
  fi
##########################################

# 이 명령어는 인터넷에서 긁어온 것을 약간 수정한 것입니다. 앞의 thread가 끝났는지 확인해서 끝났으면 바로 바로 다음 chromosome으로 넘어가도록 하려는 오더입니다.
# 중간에 sleep 10 은 10초마다 thread가 끝났는지 확인하고 만약에 앞의 thread가 멈추었다면 다음으로 넘어가게 됩니다.
# 그리고 그러한 명령이 total_chr (우리는 지금 22개의 chromosome을 설정했습니다)만큼 실행되도록 합니다. 그래서 do / done 위치를 잘 보세요. 하나의 done은 아래쪽에 있습니다.

# 그 뒤에 actual process 부분의 명령어를 & (backgroun)로 실행시키면서 thread 만큼의 process를 실행하고, 앞의 것이 끝나면 다음 것이 시행되도록 만들었습니다. 그래서 동시에 16개의 chromosome의 색인을 만들게 됩니다.

 

# Actual process of threads
# 실제로 시행되는 명령어입니다.

/usr/bin/tabix imputed/${FILE_NAME}_chr${CHR}_1K_imputed_beagle_b37_EAS.vcf.gz &
echo "Chromosome ${CHR} is being processed using Tabix."

# tabix라는 프로그램을 /usr/bin에 설치해 두었습니다. 그래서 /usr/bin/tabix로 절대 경로로명령어를 호출합니다.
# imputed/ 라는 폴더 내에, 파일들을 순서대로 indexing을 하게 됩니다. & 명령어는 background로 실행되도록 합니다. 그러면서 지금 작업중이라는 메세지를 내보내게 됩니다.

CHR=$(($CHR+1))

# CHR 1번부터 차례대로 다음으로 진행되도록 하기 위함입니다. 1을 증가시킨 뒤, 위에서 다시 total_chr와 CHR을 비교하게 됩니다.

###########################################
  done
    # 위에서 하나의 done은 아래에 있다고 했었죠?

 

lastPIDs=`jobs -l|awk '{print $2}'`
wait $lastPIDs

echo "";echo "Indexing process have been done."

# 가장 마지막 PID가 끝날 때까지 기다렸다가, 가장 마지막이 끝나면 이제 indexing이 끝났다고 알려주기 위함입니다.

 

2. Dosage filter

이제 dosage filtering을 할 겁니다. 우선 imputation된 파일을 한 번 확인해 보겠습니다. imputation을 하게 되면. vcf.gz 파일이 생성됩니다. (gz은 압축입니다.)

 

less -S Kchip_RSID_chr1_1K_imputed_beagle_b37_EAS.vcf.gz

이렇게 less로 파일을 볼 때, -S 옵션을 주면 열을 맞춰 보여줍니다.

 

##fileformat=VCFv4.2
##filedate=20231116
##source="beagle.05May22.33a.jar"
##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated ALT Allele Frequencies">
##INFO=<ID=DR2,Number=A,Type=Float,Description="Dosage R-Squared: estimated squared correlation between estim>
##INFO=<ID=IMP,Number=0,Type=Flag,Description="Imputed marker">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DS,Number=A,Type=Float,Description="estimated ALT dose [P(RA) + 2*P(AA)]">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NIH22A4207295_NIH22A4207295    
1       10177   rs367896724     A       AC      .       PASS    DR2=0.25;AF=0.3149;IMP  GT:DS   0|0:0.38     
1       10235   rs540431307     T       TA      .       PASS    DR2=0.00;AF=0.0000;IMP  GT:DS   0|0:0   0|0:0
1       10352   rs555500075     T       TA      .       PASS    DR2=0.22;AF=0.4178;IMP  GT:DS   1|0:1.09     
1       10616   rs376342519     CCGCCGTTGCAAAGGCGCGCCG  C       .       PASS    DR2=0.54;AF=0.9941;IMP

 

CHROM은 1인 것을 알 수 있고, POS는 position인데 10177, 10235 등등의 정보가 있습니다. REF와 ALT에 대한 정보도 들어가 있지요. 그러다 INFO 파트를 보시면 DR2=0.25;AF=0.3149 등등의 정보가 있는 것이 보일 겁니다. 우리는 필요한 것들만 남겨둘 것입니다.

 

DR은 imputation accuracy를 나타내는데, 0.8이나 0.9 이상이 필요하다고 보통 이야기 하는 편입니다. AF는 말 그대로 allele frequency를 의미합니다. 저는 imputation accurace > 0.9, 그리고 allele frequency > 0.1로 정의하고자 합니다. 코드는 아래와 같습니다.

 

#!/bin/bash

# name of file
FILE_NAME=Kchip_RSID

## baseline setting
chr=1
total_chr=22
thread=16

mkdir filtered_merged

##########################################
# for in. thread control
for process in $(seq 1 $total_chr)
    do j_count=`jobs -l|wc|awk '{print $1}'`
    if [[ $j_count -ge $thread ]];then
        until [[ $j_count -lt $thread ]]
            do j_count=`jobs -l|wc|awk '{print $1}'`
            sleep 10
            done
    fi
##########################################

# Actual process of threads
/home/samelsims/miniconda3/envs/bioinfo/bin/bcftools \
view -i 'DR2>0.9 & AF>0.01' imputed/${FILE_NAME}_chr${chr}_1K_imputed_beagle_b37_EAS.vcf.gz \
-Oz -o filtered_merged/${FILE_NAME}_chr${chr}_1K_imputed_filtered_beagle_b37_EAS.vcf.gz &

echo "Chromosome ${chr} is being filtered by prespecified dosage (R2)"
chr=$(($chr+1))

###########################################
    done

lastPIDs=`jobs -l|awk '{print $2}'`
wait $lastPIDs
echo "";echo "All filtering processes with cutoff value for R2 have been done."

 

위에서 설명드린 것과 나머지는 비슷 합니다. 다만 filtered_merged라는 디렉토리에 저장하기 위하여 만들었다는 점과, 실제 작동하는 코드만 좀 다를 뿐입니다.

 

# Actual process of threads
/home/samelsims/miniconda3/envs/bioinfo/bin/bcftools \
view -i 'DR2>0.9 & AF>0.01' imputed/${FILE_NAME}_chr${chr}_1K_imputed_beagle_b37_EAS.vcf.gz \
-Oz -o filtered_merged/${FILE_NAME}_chr${chr}_1K_imputed_filtered_beagle_b37_EAS.vcf.gz &

# 우선 코드 길이가 길어서 '\'로 중간에 끊어주었습니다.

# bcftools를 사용하였고, view를 통해 보면서 -i 즉, INFO column에 있는 내용들 중

# DR2>0.9, 그리고(&) AF>0.01 (MAF>1%)인 줄들만 추려내는 것입니다. DR2는 imputation accuracy라고 말씀드렸죠.

# 그리고 그 내용을, filtered라는 파일 이름을 중간에 하나 더 붙이면서 저장을 합니다.

# 그리고 위에서와 마찬가지로, & (backgroun)로 분석을 시행해서, 동시에 여러개 chromosome의 dosage filtering을 시행하게 됩니다.

 

 

이렇게 dosage filtering을 마쳤습니다. 이제 다음으로 duplication을 제거하고 하나의 plink 파일로 merge만 하면 imputation이 끝나게 됩니다.