リードのアラインメント

$ zless refseq/GCF_002575655.2.fasta.gz # 非常にサイズが大きいので、絶対に zless で開くこと

$ bwa

Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1198-dirty
Contact: Heng Li <hli@ds.dfci.harvard.edu>

Usage:   bwa <command> [options]

Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries (DEPRECATED)

         shm           manage indices in shared memory
         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ

Note: To use BWA, you need to first index the genome with `bwa index'.
      There are three alignment algorithms in BWA: `mem', `bwasw', and
      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
      first. Please `man ./bwa.1' for the manual.

$ bwa index refseq/GCF_002575655.2.fasta.gz # このコマンドは 1 日以上かかるので事前に実行済み

#!/bin/bash
 
CLEAN_FQ_DIR="clean"
BAM_DIR="bam"
 
if [ ! -d "$BAM_DIR" ]; then
        mkdir "$BAM_DIR"
else
        rm "$BAM_DIR"/*
fi
 
for ((i=1; i<=10; i++))
do
        bwa mem -t 16 -R "@RG\tID:sample${i}\tSM:sample${i}\tPL:ILLUMINA" \
                refseq/GCF_002575655.2_Aet_v5.0_genomic.fna \
                ${CLEAN_FQ_DIR}/sample${i}_R1_paired.fastq.gz \
                ${CLEAN_FQ_DIR}/sample${i}_R2_paired.fastq.gz \
                | samtools view -bF4 > ${BAM_DIR}/sample${i}.bam
done	

$ samtools view -h bam/sample1.bam | less

$ samtools view -H bam/sample1.bam

#!/bin/bash
 
BAM_DIR="bam"
CHROMOSOME_IDS=("NC_053035.2" "NC_053036.2" "NC_053037.2" "NC_053038.2" "NC_053039.2" "NC_053040.2" "NC_053041.2") # 染色体 ID
 
cd $BAM_DIR
 
for ((i=1; i<=10; i++))
do
        # alignment をフィルターする
        ## BAM ファイルをソートする
        samtools sort -@ 2 -o sample${i}_sorted.bam sample${i}.bam
        ## BAM ファイルに CSI インデックスをつける
        samtools index -@ 2 -c -o sample${i}_sorted.bam.csi sample${i}_sorted.bam
        ## BAM ファイルを filtering expression に従ってフィルタし、特定の染色体のリードのみを抽出し、ソートする
        samtools view -e '((pnext + 150) - pos) >= 70 && ((pnext + 150) - pos) <= 600' \
                -h sample${i}_sorted.bam $(IFS=" "; echo "${CHROMOSOME_IDS[*]}") | \
                samtools sort -@ 2 -o sample${i}_filtered.bam -
 
        # PCR duplicates にマークをつける
        # リード名でソート → fixmate で MC / ms タグをつける → 物理位置でソート → PCR duplicates をマークする
        samtools sort -@ 2 -n -o - sample${i}_filtered.bam | \
                samtools fixmate -m - - | \
                samtools sort - | \
                samtools markdup -@ 2 -r - sample${i}_markdup.bam
 
        rm sample${i}_sorted.bam sample${i}_sorted.bam.csi sample${i}_filtered.bam
done
 
cd -

$ ./filter-read.sh