[유전체] samtools, bedtools 이용해서 특정 부분 mapping된 fastq만들기

우선 먼저 mapping된 read만을 가져오기 위해서 아래 포스팅을 참조하세요~!

 

https://mopipe.tistory.com/150

[유전체] samtools이용해서 mapped read만 가져오기

 

 

 

JPT 샘플을 이용하여 분석해보겠습니다.

 

Target region은 chr6번을 이용하여 분석을 해보겠습니다.

 

bwa-0.7.17 mem chr6.fa.gz SRR794684_R1.fastq.gz SRR794684_R2.fastq.gz -o Target_chr6_SRR794684.sam

 

sam file을 읽어보면 다음과 같이 read별로 붙어있는 부분을 확인할 수가 있습니다.  (*표시는 붙지 못한 read)

 

여기서 samtools를 이용해서 mapped read만 가져오기 위해서는 다음과 같은 방법을 이용하면 됩니다.

 

samtools view -b -F 4 Target_chr6_SRR794684.sam >Target_chr6_SRR794684_mapped_1.sam

 

그러나 이러한 방법의 경우 mapped read 만 뽑히는 것이 아니라 단비 unmapped read를 제거하는 정도가 됩니다.

 

그래서 unmapped read만 제거된 파일을 이용하여 fastq파일을 얻으면 다음과 같은 에러가 발생이 될겁니다. 

 

pair가 맞지않은 read가 나와서 error가 발생

이것을 해결하기 위해서 mapped read + pair된 read만을 가져오기 위해서 다음과 같이 변경합니다.

 

samtools view -f 0x03 Target_chr6_SRR794684.sam -o Target_chr6_SRR794684_mapped.bam

#같은 방법으로는 다음과 같이 진행하면 됩니다.

samtools view -f 3 Target_chr6_SRR794684.sam -o Target_chr6_SRR794684_mapped_2.bam

 

위 두개의 결과값이 같은 것을 확인할 수가 있습니다.

 

이렇게 되면 mapped + pair된 read를 얻을수가 있습니다.

 

확인하시려면 다음사이트를 참고 하시면 됩니다.

 

https://broadinstitute.github.io/picard/explain-flags.html

Explain SAM Flags

 

 

이렇게 해서 mapped pair read를 가져온 bam file 확인하였고

 

bam 파일을 fastq로 만들기 위해서는 다음과 같은 tool을 이용하고자 합니다.

 

https://github.com/arq5x/bedtools2

GitHub - arq5x/bedtools2: bedtools - the swiss army knife for genome arithmetic

 

bedtool2에서는 bamToFastq라는 기능이 있어서 다음과 같이 진행해주면 됩니다.

 

bamToFastq -i Target_chr6_SRR794684_mapped.bam -fq Target_chr6_SRR794684_mapped_1.fq -fq2 Target_chr6_SRR794684_mapped_2.fq

 

그리고 나온 fastq 파일을 이용해서 분석을 하면 됩니다!

 

---