Differences

This shows you the differences between two versions of the page.

--- bioinfo:참조_서열에_대한_매핑_reference_mapping_및_시각화 [2023/06/21 14:47] – [SAM/BAM에서 unmapped read 추출하기] hyjeong
+++ bioinfo:참조_서열에_대한_매핑_reference_mapping_및_시각화 [2024/07/05 10:04] (current) – [Mapping의 실제] hyjeong
@@ Line 6: / Line 6: @@
 ===== Mapping의 실제 =====
-샘플로 사용할 Illumina sequencing read는 [[https://www.ncbi.nlm.nih.gov/sra/SRX1021861|SRX1021661]]의 run 결과 중 2백만 read(2 x 101 nt cycle)이며, read mapper로는 bowtie2를 사용한다. Bowtie2에는 paired read의 간격이나 방향 설정 등 입출력물의 형식 및 정렬 조건에 대한 다양한 옵션을 줄 수 있는데, 여기에서는 일루미나 paired end read에 해당하는 기본 조건(forward-reverse, no distance constraint)으로 하였다. Alignment의 기본 조건은 --end-to-end인데, 트리밍을 하지 않은 read라면 --local이 더 나을 수도 있다. 파라미터의 조합에 따라서 speed/sensitivity/accuracy이 달라지므로 자세한 사항은 [[https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml|bowtie2 매뉴얼]] 혹은 [[https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#getting-started-with-bowtie-2-lambda-phage-example|'Getting started' 섹션]]을 참조하라.
+샘플로 사용할 Illumina sequencing read는 //E. coli// BL21(TaKaRa)에서 유래한 [[https://www.ncbi.nlm.nih.gov/sra/SRX1021861|SRX1021661]]의 단일 run 결과(SRR2014531)로부터 2백만 read(2 x 101 nt cycle)만을 추출한 것이며, read mapper로는 bowtie2를 사용한다. Bowtie2에는 paired read의 간격이나 방향 설정 등 입출력물의 형식 및 정렬 조건에 대한 다양한 옵션을 줄 수 있는데, 여기에서는 일루미나 paired end read에 해당하는 기본 조건(forward-reverse, no distance constraint)으로 하였다. Alignment의 기본 조건은 --end-to-end인데, 트리밍을 하지 않은 read라면 --local이 더 나을 수도 있다. 파라미터의 조합에 따라서 speed/sensitivity/accuracy이 달라지므로 자세한 사항은 [[https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml|bowtie2 매뉴얼]] 혹은 [[https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#getting-started-with-bowtie-2-lambda-phage-example|'Getting started' 섹션]]을 참조하라.
   # 실습용 raw data의 설명은 https://www.ncbi.nlm.nih.gov/sra/SRR2014531 참조
@@ Line 26: / Line 26: @@
   $ samtools view -b -S -o BL21.bam BL21.sam
-SAM 파일의 대부분을 구성하는 read alignment 필드에서 두 번째 필드(flag, 0-255 범위의 숫자)는 read의 mapping 상태를 비트 단위로 나타낸다. 특정 flag를 만족하는 alignment를 추출하려면 samtools view -f INT를, 만족하지 않는 것을 추출하려면 samtools -view -F INT 명령을 사용하면 된다. Reference에 반복 서열이 존재하면 하나의 read가 여러 곳에 붙어서 각각 별도의 alignment를 생성하게 되므로(primary vs. secondary alignment), mapped read의 수치와 관련된 작업을 하려면 주의가 필요하다. 뿐만 아니라 chimeric alignment인 경우에도 하나의 read가 representative 및 supplementary alignment로 여러 차례 나타날 수 있다. SAM flags를 해독하려면 ‘samtools flags 4’와 같이 입력하거나 Decoding SAM flags 웹사이트를 이용하라. 예를 들어 0x04는 segment unmmaped를 의미한다. ‘0x’는 16진수임을 명시하기 위한 것이다(예: 10 = 0xa). [[https://gist.github.com/jvhaarst/9643139|getinsersize.py]]는 매핑 데이터를 분석하여 라이브러리 크기에 대한 정보를 제공한다.
+SAM 파일의 대부분을 구성하는 read alignment 필드에서 두 번째 필드(flag, 0-255 범위의 숫자)는 read의 mapping 상태를 비트 단위로 나타낸다. 특정 flag를 만족하는 alignment를 추출하려면 samtools view -f INT를, 만족하지 않는 것을 추출하려면 samtools -view -F INT 명령을 사용하면 된다. Reference에 반복 서열이 존재하면 하나의 read가 여러 곳에 붙어서 각각 별도의 alignment를 생성하게 되므로(primary vs. secondary alignment), mapped read의 수치와 관련된 작업을 하려면 주의가 필요하다. 뿐만 아니라 chimeric alignment인 경우에도 하나의 read가 representative 및 supplementary alignment로 여러 차례 나타날 수 있다. SAM flags를 해독하려면 ‘samtools flags 4’와 같이 입력하거나 Decoding SAM flags 웹사이트를 이용하라. 예를 들어 0x04는 segment unmmaped를 의미한다. ‘0x’는 16진수임을 명시하기 위한 것이다(예: 10 = 0xa). [[https://gist.github.com/jvhaarst/9643139|getinsersize.py]]는 매핑 데이터를 분석하여 라이브러리 크기에 대한 정보를 제공한다. 이 스크립트에 대한 상세한 설명은 [[https://allaboutbioinfo.blogspot.com/2012/04/estimating-paired-end-read-insert.html|Estimating paired-end read insert length from SAM/BAM files]]를 참조하라.
   $ samtools flagstat BL21.bam
@@ Line 126: / Line 126: @@
 자료에 따라서는 -f 12 -F 256으로 설정하기도 한다. 256(0x100)은 not primary alignment라는 뜻이므로 -F 256은 이를 부정하여 primary alignment가 된다. 그러나 unmapped read에게는 이것이 의미가 없고, 실제로 테스트한 결과 –f 12 단독으로 실행하거나 –f 12 –F 256 복합으로 실행한 경우에 차이는 없었다.
 ===== SAM/BAM에서 mapped read를 pair 형태로 추출하기 =====
+Read unmapped에 대한 SAM flag이 4(samtools view -f 0x04)이므로 이 조건에 반대되는 모든 read(-F 0x04)가 mapped read에 해당한다. 그러나 properly mapped reads, 즉 mapping된 mate의 거리와 간격이 라이브러리의 크기의 평균적 분포를 만족하는 것을 고르려면 -f 0x03(0x01 for read paired; 0x02 for read mapped in proper pair; paired read의 경우 두 flag는 동시에 쓰여야 함)을 적용해야 한다. 후속 분석 작업의 목적에 따라서는 mate 중 어느 하나만 mapping이 되었어도 read pair를 전부 추출하고 싶은 경우도 있을 것이다. SAM 및 BAM 필터링의 사례에 대해서는 [[https://gist.github.com/davfre/8596159|David Fredman의 웹페이지]]를 참조하라.
+^ **SAM flag**         ^ ** Meaning**                                                                          ^
+| [1] -f 0x04          | Unmapped                                                                              |
+| [2] -F 0x04          | Mapped                                                                                |
+| [3] -F 12            | Read and mate mapped                                                                  |
+| [4] -f 0x03          | Properly mapped \\ (-f 99/147 및 -f 83/163으로 출력되는 모든 reads)                   |
+| [5] -f 0x04 -F 0x08  | Reads that did not map, but whose mates mapped \\ (-f 4 -F 264로 표현한 자료도 있음)  |
+| [6] -f 0x08 -F 0x04  | Reads that map, but whose mates not mapped \\ (-f 8 -F 260으로 표현한 자료도 있음)    |
+여러 조건을 이용하여 복수의 BAM 파일을 추출하였다면 samtools merge 명령을 이용하여 하나로 합친 뒤 samtools bam2fq 명령을 이용하여 fastq 파일로 추출하면 된다. 위의 표에서 보인 사례에서 [5]와 [6]는 paired read까지 만들어 놓으면 같은 결과가 되므로 어떤 flag을 쓰든 관계가 없다. 다음의 사례는 BAM 파일을 탐색하여 mate 중 어느 한쪽만 mapping이 되었더라도 read pair를 추출하여 한 쌍의 fastq 파일로 출력하는 방법을 보이고 있다.
+  $ samtools view -b -F 12 BL21.bam > mapped_paired.bam
+  $ samtools view -b -f 4 -F 8 BL21.bam > unmap_map.bam
+  $ samtools view -b -f 8 -F 4 BL21.bam > map_unmap.bam
+  $ samtools merge at_least_one_mapped.bam mapped_paired.bam unmap_map.bam map_unmap.bam mapped_paired.bam
+  # Read name을 기준으로 정렬을 해야(-n) 나중에 read pair를 별도의 파일로 잘 분리해 낼 수 있다.
+  $ samtools sort -n at_least_one_mapped.bam -o at_least_one_mapped_sorted.bam
+  # -N 옵션은 read name에 /1 또는 /2를 추가한다.
+  $ samtools bam2fq –N -1 paired_1.fq -2 paired_2.fq at_least_one_mapped_sorted.bam
+추출할 read의 목록을 먼저 만든 다음에 원본 fastq 파일로부터 추출해 내는 방법도 가능하다. 이때에는 seqtk가 필요하다.
+  $ samtools view -F 12 BL21.bam | cut -f1 | sort | uniq > list_both
+  # 다음 명령어에서는 -f 4 -F 8로 실행해도 무방하다.
+  $ samtools view -f 8 -F 4 BL21.bam | cut -f1 | sort | uniq > list_single
+  $ cat list_both list_single > list_all
+  $ cp list_all list_1
+  $ cp list_all list_2
+  $ sed -i 's/$/\/1/' list_1
+  $ seqtk subseq reads_1.fq list_1 > pairedNew_1.fq
+  $ sed -i 's/$/\/2/' list_2
+  $ seqtk subseq reads_2.fq list_2 > pairedNew_2.fq
 ===== SRA data를 다운로드하는 방법(상세) =====
+본 장의 시작 부분에서 fastq-dump를 이용한 SRA data 다운로드 방법을 간략하게 설명하였다. 만약 SRA와 연계된 메타데이터 파일이 필요하거나 웹브라우저 환경의 Run Selector를 이용해야 하는 경우, 또는 아마존 웹 서비스(Amazon Web Service, AWS)를 통한 다운로드가 필요하다면 NCBI의 공식 문서인 [[https://www.ncbi.nlm.nih.gov/sra/docs/sradownload/|Download SRA sequences from Entrez search results]]를 참조하라.
+Single run에 대한 데이터를 열람하여 다운로드하려면 우선 [[https://trace.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&display=metadata|Run Browser]]에 하나의 Run accession(예: SRR2014531)을 입력하여 데이터에 접근한 다음, Reads 탭을 클릭한다. Filter 창에 조건을 입력하거나 혹은 빈 상태로 Filtered Download를 클릭한 뒤, 이어서 나타나는 창에서 Download Format을 설정한 뒤 Download를 클릭하면 된다. HTTP를 통한 다운로드이므로 한번에 하나의 데이터만 받을 수 있다.
+복수의 SRA Experiment(예: https://www.ncbi.nlm.nih.gov/sra/SRP116709)로 이루어진 데이터의 경우 원하는 Experiment에만 체크를 하여 오른쪽 위의 ‘Send to’를 클릭하여 나타나는 Run Selector로 데이터를 보낸다. 모든 Experiment를 선택하려면 체크를 하지 않은 채로 진행하면 된다. SRA Run Selector에 나타나는 Select 패널에서 Accession List를 클릭하면 SRR_Acc_List.txt 파일을 다운로드할 수 있다. 이 파일에 포함된 Run Accession을 사용하여 .sra 파일을 일괄적으로 다운로드한 뒤(SRA-Toolkit의 prefetch 사용) fastq-dump를 이용하여 fastq로 최종 전환하는 방법은 [[https://github.com/LangilleLab/microbiome_helper/wiki/Download_From_SRA|Download_From_SRA]] 문서에 설명이 되어있다. fastq-dump는 (1) NCBI에서 직접 데이터 파일(.sra)을 정해진 위치에 다운로드한 뒤 fastq로 전환하는 것과 (2) 미리 받아놓은 .sra 파일에서 fastq 파일을 추출하는 것 전부 가능하다. 다음의 예제에서는 GNU parallel을 사용하여 SRR_Acc_List.txt 파일에 수록된 Run accession을 명령어로 공급하는 기법도 쓰였다. –j 2 옵션을 주면 동시에 두 개의 커맨드를 실행하게 된다. Prefech는 ~/ncbi/public/sra/에 모든 .sra 파일을 저장하며, fastq-dump는 이로부터 fastq 파일을 추출하여 현재 위치의 fastq 디렉토리에 이를 저장한다.
+  $ PATH=/usr/local/apps/sratoolkit.2.9.6-centos_linux64/bin/:$PATH
+  $ parallel -j 1 prefetch {} ::: $(cat SRR_Acc_List.txt)
+  $ parallel -j 1 fastq-dump --skip-technical -F --split-files -O fastq {} ::: $(cat SRR_Acc_List.txt)
+SRA의 Run Browser에서 Data access 탭을 선택한 다음 URL을 클릭하여 직접 다운로드를 할 수도 있다. 아래 그림처럼 NCBI Location 오른쪽의 [[https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR8981517|URL]]을 클릭하면 SRR8981517라는 파일을 받게 된다.
+{{  :bioinfo:sra_data_access.png?400  |}}
+이를 ~/ncbi/public/sra/SRR8981517.sra라는 이름으로 저장한 뒤 다음 명령어를 써서 fastq-dump를 실행한다.
+  $ fastq-dump --split-files SRR8981517
+  Read 1431064 spots for SRR8981517
+  Written 1431064 spots for SRR8981517
+사실은 앞서 설명했듯이 ‘fastq-dump %%--%%split-files SRR8981517’라고만 입력하여 실행을 해도 .sra 파일 다운로드와 fastq 추출 작업이 자동으로 연이어서 진행된다.
 ===== 참고 자료 =====