[pysam]コンティグ末端周辺に貼りついたリードを抽出する

# コンティグ末端周辺に貼りついたリードを抽出したい

宿主ゲノム中に組み込まれたウイルスゲノムにONTリードを貼り付けた。
![igv_snapshot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/8821adae-6198-0cdc-8823-72c5bc0be085.png)
ウイルスゲノムの両末端は宿主ゲノム内の(AT/TA)nマイクロサテライトに接続している。
5'-end:
![igv_snapshot_5end.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/ffdb3cc5-9629-5046-72a7-eab0cb154f8e.png)
3'-end:
![igv_snapshot_3end.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/aa9978a2-0cb4-7c82-dd28-27e60560db8d.png)

Pythonライブラリの[pysam](https://pysam.readthedocs.io/en/latest/index.html) を使うことで、コンティグ末端に貼りついたリードから、コンティグ末端の前後数百塩基を抽出することができた。
# get_flanking_reads.py
https://github.com/satoshikawato/bio_small_scripts/blob/main/get_flanking_reads.py
## やること
コンティグ5'/3'末端から下流/上流 (window) 塩基分の領域にマップされたリードを抽出し、当該領域にマップされた配列+flank塩基分のsoft-clipped配列を出力する。
FASTAまたはFASTQ形式の出力ファイル (いずれも未圧縮) を作成するか、未指定の場合は標準出力に表示する。
## Requirements
[pysam](https://pysam.readthedocs.io/en/latest/index.html)
## Usage Example
```sh
 ./get_flanking_reads.py
usage: get_flanking_reads.py [-h] --input FILE [--output FILE] -r REF [-w WINDOW] [-m MIN] [-f FLANK] [-p PRIME] [-q]

Extract flanking reads from a BAM file

optional arguments:
  -h, --help            show this help message and exit
  --input FILE, -i FILE, --in FILE
                        Input BAM file
  --output FILE, -o FILE, --out FILE, --output FILE
                        output txt file
  -r REF, --ref REF     reference entry name
  -w WINDOW, --window WINDOW
                        window (defalt:200)
  -m MIN, --min MIN     minimum outut read lenth threshold (defalt:200)
  -f FLANK, --flank FLANK
                        flanking bases (defalt:100)
  -p PRIME, --prime PRIME
                        5'/3'-end
  -q, --fastq           output fastq
```
```sh
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR394/DRR394942/DRR394942_1.fastq.gz # Download ONT reads
seqkit seq -m 5000 DRR394942_1.fastq.gz -o DRR394942.up5k.fq.gz # Remove <5kb reads
# Align length-filtered ONT reads against the assembly, select only well-mapped reads, then realign the selected reads
NUMTH=8
LIB=LC738884
N=5
LONGR=DRR394942.up5k.fq.gz
(minimap2 --MD -ax map-ont -t $NUMTH ${LIB}.fa ${LONGR}| samtools view -Sb -q 60 -  |samtools sort -@ $NUMTH -O BAM -o ${LIB}.up${N}k.tmp.bam)2>${LIB}.up${N}k.tmp.minimap2.log
  samtools index ${LIB}.up${N}k.tmp.bam
  samtools view -h ${LIB}.up${N}k.tmp.bam| paftools.js sam2paf -p - >${LIB}.up${N}k.tmp.paf
awk '{FS="\t"; OFS="\t"};$10/$11 >= 0.8 && $11 >= 5000 {print $1}' ${LIB}.up${N}k.tmp.paf |sort -u > ${LIB}.up${N}k.tmp.mapped.list
samtools fastq ${LIB}.up${N}k.tmp.bam |gzip >${LIB}.up${N}k.tmp.fq.gz
seqkit grep -f  ${LIB}.up${N}k.tmp.mapped.list ${LIB}.up${N}k.tmp.fq.gz | gzip >${LIB}.up${N}k.mapped.fq.gz
(minimap2 --MD -ax map-ont -t $NUMTH ${LIB}.fa ${LIB}.up${N}k.mapped.fq.gz| samtools view -Sb -q 60 -  |samtools sort -@ $NUMTH -O BAM -o ${LIB}.up${N}k.bam)2>${LIB}.up${N}k.minimap2.log
  samtools index ${LIB}.up${N}k.bam

./get_flanking_reads.py -i LC738884.up5k.bam -r LC738884.1 -o LC738884.up5k.5end.fa -w 100 -m 200 -f 100 -p 5
./get_flanking_reads.py -i LC738884.up5k.bam -r LC738884.1 -o LC738884.up5k.3end.fa -w 100 -m 200 -f 100 -p 3
```

```LC738884.up5k.5end.fa
>DRR394942.67641:3365-3565
TATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAATATATATATATATATATGTTGTAGCTCTAATAAGGATAATGACTCATTAACGCTAATAGGTAACGCATTTTATTAGGCAGGGAAAACTGTCCATTAACGTTAATGGACAGTTTTGCC
>DRR394942.72653:3794-3994
ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATTTATATATATATATATATATATATATATATATATGTTGTGCCTAATAAGGATAATGACTCATTAACGCTAATAGGTACTCCATTTTTATTAGGCAGGAAAACTGTCCATTAACGTTAATGGACAGTTTTGCCTG
>DRR394942.150247:911-1111
ATATATATATATATATATATACATACATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTTGTGCCTAATAAGGATAATGACTCATTAACGCTAATAGGTAACGCATTTTTTACAGAAAACTGTCCATTAACGTTAATGGACAGTTTTGCCTGCATAA
>DRR394942.172995:8555-8755
ATATGCATATATATAAATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTTGTGCCTAATAAGGATAATGACTCATTAACGCTAATAGGTAACGCATTTTTTATTAGGCAGGGAAAACGTCCATTAACGTTAATGGACAGTTTTGCCT
>DRR394942.183038:960-1160
ATATATATATATATATATAAATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTTGTGCTTAATAAGGATAAAGACTCATTAACGCTAATAGGTAACGCATTTTTTATTAGGCAGGGAAAACTGTCCATTAACGCATGGACAGTTTTGCCTG
>DRR394942.187484:8678-8878
ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTTGTGCCTAATAAGGATAATGACTCATTAACGCTAATAGGTAACGCATTTTTTATTAGGCAGGAAAAAACTGTCCATTAACGTTAATGGACAGTTTTGC
>DRR394942.188789:1576-1776
ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTTGTGCCTAATAAGGATAATGACTCATTAACGCTAATAGGTAACGCATTTTTTATTAGCAGGGAAAACTGTCCATTAACGTTAATGGACAGTTTTGCCT
>DRR394942.325:4953-5153
ATATATATATATATACATACATATATATATATATATATATATATATATATATATATGAAGCCCTGCTACACAATGACCTACTAACGTCAATAGCAACGACATTAGGCAGGGAAACCAATCCATTAACGTTAATGGCAGTTTATTCACATGAATAAGACAAGAAGAGTCCTATTGTGGTGATTCATTATATCAATGGGGTT
```
```LC738884.up5k.3end.fa
>DRR394942.163225:25582-25782
AGTGTTTCATAGTGGGCAATGTATTCCGTCTAGTATGATTAATAGGACAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTATTAGGCACAACATATATATATATATATAATATATAATATATATATATATATATATATATATATATATATATATACATATATATATATTTTTTTTTTTCTCTCTCTCTCTCT
>DRR394942.114659:21151-21351
GTGTTTCATAGTGGGCAATGTATTCCGTCTAGTATCATTAATGGACAAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTTATTAGGCACAACATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAA
>DRR394942.141448:20346-20546
GTGTTTCATAGTGGGCAATGTATTCCGTCCTAGTATGATTAATGGACAAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTATTAGGCACAACATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATTACATTTACATATATATATATATATATATATATA
>DRR394942.168575:17703-17903
AGTGTTTCATAGTGGGCAATGTATTCCGTCTAGTATGATTAATGGACAAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTATTAGGCACAACATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATTATATATATATATATATATATATATATATATATA
>DRR394942.148858:13288-13488
AGTGTTTCATAGTGGGCAATGTATTCCGTCTAGTATGATTAATGGACAAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTATTAGGCACAACATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAT
>DRR394942.182075:9927-10127
GTAAGTGTTCCTTAGTGGGCAATGTATTCCGTCTAGTATGATTAATGGACAAAGTGTCGTTACCTATTAGCGTTGAGTCATTATCCTTATTAGGCACAACACATATATATATATATTATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATATGTCTTTATATATGTATA
>DRR394942.98765:9426-9626
CAGTGCATCATAGTGGGCAATGTATTCAATCTAGTATGATTAATGGACAAGTGTCGTTACCTATTAGCGTTAATGAGTCATTATCCTTATTAGGCACAACATATATATATATATATATATATATATATATATATATATATATATATATATATTGTATATCTATATGTGTATATATATATATATATATATATATATATATA
```
## Code
pysamではFastqProxyオブジェクトを直接インスタンス化できないので新たに```ReadObject```クラスを作った。

```python:get_flanking_reads.py
#!/usr/bin/env python
# coding: utf-8

import sys,os
import pysam
import argparse
from contextlib import redirect_stdout

def _get_args():
    parser = argparse.ArgumentParser(description='Extract flanking reads from a BAM file')
    parser.add_argument("--input","-i", "--in",metavar="FILE",help="Input BAM file",required=True)
    parser.add_argument("--output","-o","--out", "--output", metavar="FILE", help="output txt file")
    parser.add_argument("-r","--ref",type=str,help="reference entry name",required=True)
    parser.add_argument("-w","--window",type=int,help="window (defalt:100)",default=100)
    parser.add_argument("-m","--min",type=int,help="minimum outut read lenth threshold (defalt:200)",default=200)
    parser.add_argument("-f","--flank",type=int,help="flanking bases (defalt:100)",default=100)
    parser.add_argument("-p","--prime",type=int,help="5'/3'-end",default=5)
    parser.add_argument("-q","--fastq",help="output fastq",action='store_true')
    if len(sys.argv) == 1:
        parser.print_help(sys.stderr)
        sys.exit(1)
    args = parser.parse_args()
    return args

class ReadObject:
    def __init__(self, read_name, read_sequence, read_quality):
        self.name = read_name
        self.sequence = read_sequence
        self.quality = read_quality

def calculate_read_coordinates(read, window, flank, prime):
    if prime == 5:
        read_start = read.query_alignment_start - flank
        read_end = read.query_alignment_start + window 
        read_length = len(read.query_sequence[read_start:read_end])
    else:
        read_start = read.query_alignment_end - window
        read_end = read.query_alignment_end + flank 
        read_length = len(read.query_sequence[read_start:read_end])
    return read_start, read_end, read_length
    
def create_read_object(read, read_start, read_end):
    new_read_name = "{}:{}-{}".format(read.qname, read_start, read_end)
    new_read_sequence = read.query_sequence[read_start:read_end]
    new_read_quality = str(''.join(map(lambda x: chr( x+33 ), read.query_qualities[read_start:read_end])))
    read_object = ReadObject(new_read_name, new_read_sequence, new_read_quality)
    return read_object

def print_fastq_entry(read):
    print('@{}'.format(read.name))
    print(read.sequence)
    print("+")
    print(read.quality)
    
def print_fasta_entry(read):
    print('>{}'.format(read.name))
    print(read.sequence)

def main():
    args = _get_args()
    in_bam = args.input
    out_file = args.output
    flank= args.flank
    window = args.window
    read_length_threshold = args.min
    prime=args.prime
    ref=args.ref
    fastq = args.fastq
    bamfile = pysam.AlignmentFile(in_bam, "rb")
    reference_length = bamfile.get_reference_length(ref)
    if prime == 5:
        start = 1
        end = window
    else: 
        start = reference_length - window
        end = reference_length
    reads = list()
    for read in bamfile.fetch(ref, start, end):
        read_start, read_end, read_length = calculate_read_coordinates(read, window, flank, prime)
        if read_length < read_length_threshold:
            continue
        read_object = create_read_object(read, read_start, read_end)
        reads.append(read_object)
    if out_file is not None:
        with open(out_file, 'w') as f:
            with redirect_stdout(f):
                for read in reads:
                    if fastq == True: 
                        print_fastq_entry(read)
                    else: 
                        print_fasta_entry(read)
    else:
        for read in reads:
            if fastq == True: 
                print_fastq_entry(read)
            else: 
                print_fasta_entry(read)

if __name__ == "__main__":
    main()
```

# References
https://pysam.readthedocs.io/en/latest/index.html