[SAMtools] BAM에서 특정 region에 대한 fastq 뽑기

SAMtools를 이용하면 특정 region에 대한 fastq를 얻어 낼 수 있다. 

준비물
BAM 
특정 region에 대한 BED 

samtools view -b -L [BED] -@ 10 [BAM]| samtools fastq -1 [paired 1 fastq] -2 [paired 2 fastq] -




근데 해당 분석을 진행하면, paired end에서 BED region에 overlapping되는 부분만 얻을 수 있게 된다. 

fastq1과 fastq2의 pair가 맞지 않는 상황이 발생한다.

예) 30X로 만들어진 BAM, 유전자 한개(2 kb 정도)에 대하여 수행한 결과

Pair 1은 188 reads

Pair 2는 195 reads이고

Pair 양쪽에 있는 경우는 163개였다.

합집합은 220 이기 때문에 paired 갖춘 fastq만 뽑을경우 163 paired 만 사용하게 되므로 163 * 2 / (188 + 195) = 85% 의 데이터만 이용하게 된다.

그렇기 때문에

1) 그냥 single end로 BWA MEM을 진행한다.

2) 한쪽 fastq만 얻어진 QNAME에 대하여 bam에서 다시 fastq를 추출한다

3) intersect만 분석한다.

라는 세 가지 방법이 있게 된다.

2의 경우에는 손이 많이 가기 때문에 우선 1로 분석해보자!

3도 쉽긴 한데, 15% 데이터 손실이 있게 된다.