BWA-aln简介
- BWA（Burrow-Wheeler Aligner），是一款将DNA序列mapping到参考基因组上的软件，能够将差异度较小的序列比对到一个较大的参考基因组上。有三个比对算法：BWA-backtrack，BWA-SW和BWA-MEM。
- BWA-aln对应的就是BWA-backtrack算法：用于将100bp以内的短序列比对到参考基因组（本来是用来比对 Illumina 的序列的，reads 长度最长能到 100bp），此处就借用来将多条50bp的序列比对到基因组。
- BWA-aln对应的命令：bwa aln/samse/sampe（samse=sample single-end，sampe=sample paired-end）。
项目需求
- 将上千个50bp的染色体探针序列比对到最新版基因组
序列准备
- FASTA格式探针文件（BWA软件支持格式为fasta、fastq以及他们的压缩文件.gz作为比对序列）,可以使用AWK进行调整
- 基因组FASTA文件,可以使用
$ less genome.fa | awk '$1 ~"^>"'查看所有序列名（包括染色体），有的基因组含有未组装的Scaffold序列，可以拼作为一条未知染色体，还有的基因组含有线粒体、叶绿体基因组等。
基因组比对
- BWA安装：注意配置 .conda 文件
$ conda create -n bwa python=3.7
$ conda activate bwa
$ conda install bwa
- 建立基因组索引 → 比对到基因组
# bwa建立索引
bwa index csi.chromosome.fa
# 比对到基因组，生成二进制文件（-t 指定线程数）
bwa aln -t 2 -f ./probe1.sai ./genome.fa ./probe1.fa
# 转为sam文件，因为是一条序列的比对，使用单端测序的解读方式
bwa samse -f ./probe1.sam ./genome.fa ./probe1.sai
# .sai文件可以删除了
探针比对结果统计
-
sam文件的解读：生信：2：sam格式文件解读，博主写的很全，很厉害
- 主体部分有11个主列和1个可选列
QNAME 比对的序列名称 例如：M04650:84:000000000-B837R:1:1101:22699:1759（一条测序reads的名称）
FLAG Bwise FLAG（表明比对类型：paring，strand，mate strand等） 例如：99
RENAME 比对上的参考序列名 例如：NC_000075.6
POS 1-Based的比对上的最左边的定位 例如：124057649
MAPQ 比对质量 例如：60
CIGAR Extended CIGAR string（操作符：MIDNSHP）比对结果信息；匹配碱基数，可变剪接等 例如：87M
MRNM 相匹配的另外一条序列，比对上的参考序列名 例如：=
MPOS 1-Based leftmost Mate Position （相比于MRNM列来讲意思和POS差不多） 例如：124057667
ISIZE 插入片段长度 例如：200
SEQ 和参考序列在同一个链上比对的序列（若比对结果在负义链上，则序列是其反向重复序列，反向互补序列） 例如：ATTACTTGGCTGCT
QUAL 比对序列的质量（ASCII-33=Phred base quality）reads碱基质量值 例如：-8CCCGFCCCF7@E-
可选的列 以TAG：TYPE：VALUE的形式提供额外的信息
- 主体部分有11个主列和1个可选列
写脚本统计sam结果
$ cat stat2.sh
#! /bin/sh
#
cd $(pwd)
echo -e "probe\tchr1\tchr2\tchr3\tchr4\tchr5\tchr6\tchr7\tchr8\tchr9" > title.txt
for i in $(ls *.sam)
do
echo ${i%.*} > ${i%.*}.prob
for j in 'chr1' 'chr2' 'chr3' 'chr4' 'chr5' 'chr6' 'chr7' 'chr8' 'chr9'
do
less $i | awk '$3 == "'$j'" {print $0}' | wc -l >> ${i%.*}.prob
done
done
for i in $(ls *.prob)
do
paste -s $i > tmp.row
cat title.txt tmp.row > probedistribution.csv
cat probedistribution.csv > title.txt
done
- 统计效果
$ cat probedistribution.csv
probe chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9
chr1a 1200 0 0 0 0 0 0 0 0
chr1b 1003 0 0 0 0 0 0 0 0
chr2a 0 1292 0 0 0 0 0 0 0
chr2b 0 1685 0 0 0 0 0 0 0
chr3a 0 0 1065 0 0 0 0 0 0
chr3b 0 0 1097 0 0 0 0 0 0
chr4a 0 0 0 1330 0 0 0 0 0
chr4b 0 0 0 1047 0 0 0 0 0
chr5a 0 0 0 0 1344 0 0 0 0
chr5b 0 0 0 0 1177 0 0 0 0
chr6a 0 0 0 0 0 875 0 0 0
chr6b 0 0 0 0 0 832 0 0 0
chr7a 0 0 0 0 0 0 1106 0 0
chr7b 0 0 0 0 0 0 1110 0 0
chr8a 0 0 0 0 0 0 0 849 0
chr8b 0 0 0 0 0 0 0 1112 0
chr9a 0 0 0 0 0 0 0 0 1134
chr9b 0 0 0 0 0 0 0 0 1182
后续
- 获取bed文件后得到密度文件
- RIdeogram绘制染色体探针分布图(https://www.jianshu.com/p/9b83122a9cdf)
