NGS 데이터 분석 - CRAM 파일 핸들링

• 현재 나에게 FASTQ & CRAM & VCF 파일이 있지만, FASTQ는 sequencing에 대한 어떠한 정보가 담겨 있지 않아 uBAM 파일로 만들어 봤자 염기서열만 존재하는 깡통이다.
• 따라서 이미 Alignment가 된 CRAM 파일을 활용해야 하는데, GATK 전체 파이프라인 이해를 위하여 BAM으로 다시 변환하고자 한다.

이때, mapping 당시의 reference genome이 필요한데 이에 대한 정보가 전혀 없다.

 

[1] CRAM header 추출

samtools view -H input.cram

# UR:/mnt/ssd/MegaBOLT_scheduler/reference/hg38.fa
# @PG ID:bwa PN:bwa VN:BWA_pre2.0.1.2 CL:/mnt/ssd/MegaBOLT_scheduler/bin/bwa mem /mnt/ssd/MegaBOLT_scheduler/reference/hg38.fa -l /mnt/ssd/MegaBOLT_scheduler/tmpDir/14/1/NB4Q6634/megabolt.inputlist -t 40 -i /mnt/ssd/MegaBOLT_scheduler/lic/MegaBOLT.lic
# @PG ID:bamsormadup PN:bamsormadup CL:/mnt/ssd/MegaBOLT_scheduler/bin/SortMarkDup threads=10 tmpfile=/mnt/ssd/MegaBOLT_scheduler/tmpDir/14/1/tmp_sort_5 inputformat=sam outputformat=bam indexfilename=/mnt/ssd/MegaBOLT_scheduler/tmpDir/14/1/NB4Q6634.bwa.sortdup.bam.bai PP:bwa VN:2.0.79
# @PG ID:samtools PN:samtools PP:bamsormadup VN:1.10 CL:/mnt/ssd/MegaBOLT_scheduler/bin/samtools view -SC -T /mnt/ssd/MegaBOLT_scheduler/reference/hg38.fa -o /mnt/ssd/MegaBOLT_scheduler/tmpDir/14/1/NB4Q6634.bwa.sortdup.cram --threads 10 --write-index /mnt/ssd/MegaBOLT_scheduler/tmpDir/14/1/NB4Q6634.bwa.sortdup.bam

• Reference로 'hg38.fa' 경로를 지정하였고, bwa, bamsormadup, samtools를 사용하였다.

[2] Reference 찾기

(1) GATK recource bundle

• CRAM 파일의 header를 보면 MegaBOLT를 활용한 것으로 나타나는데, manual (pdf)을 보니 GATK resource bundle (ftp, link)에서 다운받은 것으로 명시되어 있다.

ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/ # no password

(2) NCBI

• 나와 같은 상황인 사람이 Reddit에 쓴 글(link)에 따르면, 이 사람의 샘플은 hg38 without ALT/HLA contigs or decoys (ftp)을 사용한 것으로 보인다고 한다.

Looks like hg38 was used as the reference

EDIT: now that I have my results, I think I can confirm that the other poster had it correct, and GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz is the reference sequence.

For what it's worth, I've just been told by a Nebula customer support agent that they use the hs38 analysis set without ALT/HLA contigs or decoys, in case you were specifically curious about the reference used.

 

• hg38를 사용한 것은 맞지만, 정확한 파일을 찾아야한다. NCBI ftp site에 들어가보니 (link) 다양한 버전의 파일이 존재한다.
  • GCA_000001405.15_GRCh38_full_analysis_set.fna.gz
  • GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

[3] CRAM과 Reference 비교

위에서 2가지 후보군으로 추려서 각 버전을 다운받았다.

• GATK bundle set
• NCBI no ALT/HLA

diff --side-by-side \
<(samtools view -H input.cram | cut -f 1-4) \
<(samtools dict ref.fna | cut -f 1-4)

(좌) input.cram (우) hg38 reference

• 위의 그림과 같이 hg38과 일치하는 것을 확인하였다. 대신, 아래와 같이 ALT/HLA region은 없다.

(좌) input.cram (우) ALT/HLA contigs or decoy

• 따라서 Reddit에 언급되었던 버전이 alignment에 사용된 reference으로 확인되었다.

[4] CRAM에서 BAM 변환

samtools view -b \
  --reference ref.fna \
  -o output.bam
  input.cram

• 약 5시간이 걸려 CRAM에서 BAM 파일로 변환하였다.