NGS 데이터분석 1 -예제 데이터로 분석해보기(SRR8241280~1)

데이터분석 공부도 하고 GATK pipeline도 만들겸 

public data로 분석해보기로 한다.

포스팅 1편에서는 GATK를 사용하기 전 단계인 전처리 단계를 진행한다.

 

 

SAMPLE 

내가 사용한 데이터는 project ID PRJNA505117의 두개의 sample이다.

 

Study Title은 Whole exome sequencing of Colorectal Cancer Liver Metastases이다.

 

SRR8241280~SRR8241309 의 많은 sample들이 있지만 전부 돌려보기엔 시간과 컴퓨터 자원의 한계로 인해

SRR8241280,SRR8241281 두개의 sample을 분석해본다.

 

SRR8241280: TUMOR sample

SRR8241281: NORMAL sample

 

 

본격적으로 분석을 진행하기 전에 somatic variants calling 하기위한 기본 input 준비물을 정리해본다.

---

DATA

1. reference 시퀀스 파일 : Homo_sapiens_assembly38.fasta

reference는 GATK 에서 제공하는 Google Cloud에서 다운 받았다.

GATK에서 indexing까지 해놓은 파일들까지 다운 받을 수 있다. 하지만 내가 사용하는 BWA 버전이랑 다르면 error를 만날수 있다. 

 

console.cloud.google.com/storage/browser/genomics-public-data/references/hg38/v0?pageState=(%22StorageObjectListTable%22:(%22f%22:%22%255B%255D%22))&prefix=&forceOnObjectsSortingFiltering=false

 

Google Cloud Platform

# ncbi와 같은 db에서도 쉽게 reference seq data를 받을 수 있다. 다만 unmapped contigs들이 들어있을 수 있기때문에 직접 parsing을 해주어야한다.

 

2. read 파일 : 위에서 설명한 SRR8241280_1, SRR8241280_2, SRR8241281_1, SRR8241281_2

                  파일명의 _1, _2는 paired read쌍을 의미한다. -> sample두개의 각각 한쌍 

 

3. Pannel of Normal (PoN) : 1000g_pon.hg38.vcf.gz, 1000g_pon.hg38.vcf.gz.tbi

PON은 체세포 변이 분석에 사용되는 리소스 유형이다.  recurrent technical artifacts를 제거하는것이 주 목적이기때문에 같은 플랫폼으로 만든 PON을 사용해야 한다. 1000지놈으로 만든 PON은 GATK에서 제공한다. 나는 indexing된 파일과 함께 받았다.

 자신의 PON도 만들수있다. 다만 최소 40개의 정상 샘플이 필요하다.

What all PONs have in common is that (1) they are made from normal samples (in this context, "normal" means derived from healthy tissue that is believed to not have any somatic alterations) and (2) their main purpose is to capture recurrent technical artifacts in order to improve the results of the variant calling analysis.

출처 : gatk.broadinstitute.org/hc/en-us/articles/360035890631-Panel-of-Normals-PON-

console.cloud.google.com/storage/browser/gatk-best-practices/somatic-hg38;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

Google Cloud Platform

4. germline-resource : af-only-gnomad.hg38.vcf.gz,  af-only-gnomad.hg38.vcf.gz.tbi

 

위의 PoN을 받은 곳에서 마찬가지로 indexing된 파일과 함께 받아준다.

 

data 준비 완료.

---

 

ANALYSIS

원래는 QC단계가 포함되어야 하지만 public data라 포함하지 않았다.

추후 추가 할 예정이다.

 

 

변수 설정

REF="/data/path/is/Genomics/data/ref/Homo_sapiens_assembly38.fasta"
DATA_PATH='/data/path/is/Genomics/data/'
RESULT_PATH='/data/path/is/Genomics/result/'
RESULT_DIR1='1_bwa'
N_THREADS=25

RESULT_BWA=${RESULT_PATH}${RESULT_DIR1}

PICARD="/data/path/is/Genomics/tools/picard.jar"
SAMTOOLS="/data/path/is/Genomics/tools/samtools-1.11/samtools"
ALL_STARTTIME=$(date +%s)

BWA mem

RG을 지정해줘야 추후에 분석을 진행하는데 편하다고 한다.

-R옵션을 준다.

bwa index -a bwtsw $REF
bwa mem -t $N_THREADS -R "@RG\tID:foo\tLB:bar\tPL:illumina\tPU:illumina\tSM:${1:0:10}" $REF ${1} ${2} > ${RESULT_BWA}/${1:0:10}.sam
bwa mem -t $N_THREADS -R "@RG\tID:foo\tLB:bar\tPL:illumina\tPU:illumina\tSM:${3:0:10}" $REF ${3} ${4} > ${RESULT_BWA}/${3:0:10}.sam

SAMtoBAM

$SAMTOOLS view -bhS -@${N_THREADS} $SAM > ${RESULT_BWA}/${NAME:0:-4}.bam

Sorting BAM

$SAMTOOLS sort $BAM -o ${RESULT_BWA}/${NAME:0:-4}.sorted.bam -@ ${N_THREADS}

Markduplicates

        java -jar $PICARD MarkDuplicates \
        I=$FILE \
        O=${RESULT_BWA}/${NAME:0:-4}_dedup.bam \
        METRICS_FILE=duplicates \
        REMOVE_DUPLICATES=true \
        CREATE_INDEX=true

Reference dictionary & indexing

java -Xmx1024m -jar $PICARD CreateSequenceDictionary \
R=$REF \
O=${REF:0:-6}.dict

$SAMTOOLS faidx $REF

GATK를 다뤄보기 위한 준비가 끝났다.

 

2편에서는 GATK를 다뤄 보쟈.