Reference sequence에 매핑된 read를 뽑아내는 방법은 몇 가지가 존재한다. CLC Genomics Workbench에서도 다음 사이트를 참조하여 실행하면 된다.

CLC Exporting Mapped Reads

그러나 이 방법이 특정 reference sequence(혹은 특정 region)에 mapping된 read만을 꺼낼 수 있는지, 혹은 reference set 전체에 대해서도 적용 가능한지는 자세히 확인하지 않았다.

이 페이지에서는 read pair 중 어느 하나만 reference에 mapping이 되었어도 pair를 전부 추출하여 후속 분석 작업에 활용하는 방법을 알아보고자 한다. 다음 그림을 참조하면 쉽게 이해가 될 것이다.

사용할 read mapper는 Sanger Institute의 SMALT이다.

SMALT employs a hash index of short words up to 20 nucleotides long and sampled at equidistant steps along the reference genome. For each sequencing read, potentially matching segments in the reference genome are identified from seed matches in the index and subsequently aligned with the read using dynamic programming.

SMALT v0.7.6은 bioconda 환경으로 설치해 두었다(base environment).

$ smalt index -k 14 -s 8 REF reference.fasta
$ smalt map -n 16 -f bam -o mapped.bam REF reads_1.fastq reads_2.fasta  
$ samtools stats mapped.bam | grep ^SN | cut -f 2- # mapping report 출력

-k 14 -8의 의미는 reference.fasta 파일에서 길이 14 bp의 word를 매 8번째 위치마다 샘플링한다는 뜻이다. 이렇게 인덱스 처리된 reference sequence는 “REF”라는 이름으로 mapping에 사용하면 된다.

나의 질문은 Biostars의 다음 질문과 매우 흡사하다.

1. Extract One-End Mapped Paired-End Reads From Bwa Bam/Sam Files
2. Collect Read Pairs Where At Least One Read Is Mapped
3. (SeqAnswers) SAM: Extract both pairs where at least one pair mapped

다음의 방법은 유일한 것도 아니고 가장 효율적인 것도 아닐 것이다. 그러나 원칙에 가장 충실한 방법일 것이라 자부한다. Read 목록을 참조하여 원본 read file에서 다시 서열을 뽑아내는 유틸리티로서 seqtk를 사용하였다. 이것 역시 bioconda로 재설치하였다.

$ samtools view -F 12 mapped.bam | awk 'NR%2!=0{print $1}' > list_both
$ samtools view -f 8 -F 4 mapped.bam | awk '{print $1}' > list_single
$ cat list_both list_single > list_all
$ seqtk subseq reads_1.fastq list_all > subset_1.fastq
$ seqtk subseq reads_2.fastq list_all > subset_2.fastq

첫번째 명령에서는 samtools view의 출력물 중에서 홀수 라인만을 추출하였다. 이렇게 하지 않으면 read name이 연달아 두 줄에 나타나기 때문이다(pair이므로) ⇐ 그런데 지금 회고해 보니 정확하지 않음을 알 수 있었다. 왤까? — Haeyoung Jeong 2019/08/07 14:41