我就是Super Star——基因定位之BSA

前两天有小伙伴留言介绍BSA，这里继续带大家回顾下李广伟师兄在2017年发的BSA教程，之后我还会写一篇有关的文章讲讲BSA实战

话说到了基因组时代，通过正向遗传学，找到一个基因，基本就是三个方法：GWAS、bin-map、BSA。无论何种方法，用的基本原理还是我们高中生物课本上学的：

1. 基因在染色体上直线排列，

2. 染色体会发生重组和交换。

基因组很大，几万个基因，我们很难直接找到影响某个表型的cause基因。今天要讲的就是BSA。

先给大家讲一个故事

假如我们有两个主角“金角大王”和“银角大王”。金色还是银色，只有一个基因控制。两位大王基因组上，除了我们目标基因不同造成颜色的差异，还有三个SNP，分别为SNP1、2、3（如图所示）。
注意：原来金色allele和A，G，C在一条染色体上，银色和T，C，G在同一条染色体上。

下面故事发生了：
金角大王和银角大王结婚了，然后生了一群小大王，小大王自由婚配生了一群小小大王…….一群小妖，有金色的也有银色的。有一天，我们抓了一群金色小妖，测个序，看了下每个位点上等位基因频率。发现这个金色小妖群里SNP2上G的等位基因频率很高，SNP1上A次之，SNP3上的C的最低（注意到：开始时候SNP2G，SNP1A，SNP3C都是和金色相随相伴的，即起始等位基因频率都是1）。

那么，如果我们只知道SNP1、2、3的等位基因频率，而不知道金色银色基因在哪的话，我们可以基本判断目标基因在SNP2附近，而不在SNP3附近。

为什么呢？因为越近、越连锁，关联性越强。我们选择这些极端表型的时候，把这个控制表型目标基因的邻居顺便都捞出来了，然后我们可以通过邻居找到当家的。就是这么简单粗暴的大道理。

也就是说，我们可以通过表型选择，把金角小妖混一块，银角小妖混一块，然后扫描全基因组所有位点等位基因频率差异，找到颜色基因所处的位置。

这就是BSA的基本原理：

闲话少说，最近我崇拜的大神，冷泉港的lippman在Cell上灌了一篇封面，两个主要的基因都是通过BSA进行定位，然后用同源基因方法找到的，具体信息自己看文章。文章链接：http://www./cell/fulltext/S0092-8674(17)30486-5
为了表示对大神的崇拜之情，我第一时间把文章的信息down下来，进行了重演，以下就是详细过程。
所用软件版本：

• bwa:0.7.9a-r78

• samtools:0.1.19-44428cd

• bcftools:0.1.19-44428cd

实际流程

下载数据和软件安装

根据文章 NCBI SRA的索取号 下载

wget -c ftp://ftp-trace.ncbi./sra/sra-instant/reads/ByRun/sra/SRR/SRR527/SRR5274882/SRR5274882.sra
wget -c ftp://ftp-trace.ncbi./sra/sra-instant/reads/ByRun/sra/SRR/SRR527/SRR5274881/SRR5274881.sra
wget -c ftp://ftp-trace.ncbi./sra/sra-instant/reads/ByRun/sra/SRR/SRR527/SRR5274880/SRR5274880.sra

需要安装软件 SRA tools

wget https://ftp-trace.ncbi.nlm./sra/sdk/2.8.2-1/sratoolkit.2.8.2-1-centos_linux64.tar.gz
tar zxvf sratoolkit.2.8.2-1-centos_linux64.tar.gz
##直接解压缩即可得到可执行文件
### test [ligw@node13 sratoolkit.2.8.2-1-centos_linux64]$ ./bin/fastq-dump -h

NCBI下载的SRA格式转换成fastq格式

/home/ligw/tools/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump -F --skip-technical --split-files  --defline-qual '+' --defline-seq '@$ac-$si/$ri' --split-3 SRR5274880.sra&
/home/ligw/tools/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump -F --skip-technical --split-files  --defline-qual '+' --defline-seq '@$ac-$si/$ri' --split-3 SRR5274881.sra&
/home/ligw/tools/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump -F --skip-technical --split-files  --defline-qual '+' --defline-seq '@$ac-$si/$ri' --split-3 SRR5274882.sra&

参考基因组：

 wget ftp://ftp.ensemblgenomes.org/pub/plants/release-35/fasta/solanum_lycopersicum/dna/Solanum_lycopersicum.SL2.50.dna.toplevel.fa.gz

SNP calling pipeline ：

### index  
bwa index tomato.SL2.50.fa
### mapping    
bwa mem -M -t 20 -R '@RG\tID:ejmt\tSM:ejmtpool\tPL:Illumina' tomato.SL2.50.fa SRR5274882_1.fastq SRR5274882_2.fastq > ejmt.sam
### convert sort index remove_duplicate index
samtools view -bS -o ejmt.bam ejmt.sam
samtools sort -m 2G -@ 20 ejmt.bam ejmt.sorted
samtools index ejmt.sorted.bam
### MarkDuplicates
java -Xmx20G -jar /home/ligw/tools/picard.jar MarkDuplicates I=ejmt.sorted.bam O=ejmt.sorted.redup.bam \
METRICS_FILE=ejmt.sort.metrics REMOVE_DUPLICATES=true \
ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT

samtools index ejmt.sorted.redup.bam
## vaiant calling
samtools  mpileup -DSug -C 50 -Q 20 -q 40 -f tomato.SL2.50.fa ejmt.sorted.redup.bam  | bcftools view -cvg - > ejmt.vcf

PS: 野生型池就不重复了

Mapping 基因

shell部分：

### 过滤掉indel

grep -v '##' ejmt.vcf | grep -v 'INDEL' > ejmt.snp.vcf
grep -v '##' jmt.vcf | grep -v 'INDEL' > jmt.snp.vcf
grep -v '##' wt.vcf | grep -v 'INDEL' > wt.snp.vcf

R语言代码：

e.snp <- read.table(file='ejmt.snp.vcf',head=T,as.is=T)
e.snp <- e.snp[e.snp$CHROM==1|e.snp$CHROM==2|e.snp$CHROM==3|e.snp$CHROM==4|e.snp$CHROM==5|e.snp$CHROM==6|
e.snp$CHROM==7|e.snp$CHROM==8|e.snp$CHROM==9|e.snp$CHROM==10|e.snp$CHROM==11|e.snp$CHROM==12,]

#### 利用VCF文件里，INFO一栏DP4信息，得到ref allele read和 alt allele reads#####

index <- lapply(e.snp$INFO,function(x){  
   str <- regexpr('DP4',x)[1]  end <- regexpr('MQ',x)[1]  
   DP4.str <- sub('DP4=','',substr(x,str,end-2))  
   return(DP4.str)})

ref.no <-  unlist(lapply(index,function(x){
   sum(as.numeric(unlist(strsplit(x,',')))[1:2])}))

alt.no <-  unlist(lapply(index,function(x){  
   sum(as.numeric(unlist(strsplit(x,',')))[3:4])}))

e.snp.index <- e.snp[,1:6]
e.snp.index$ref.no <- ref.no
e.snp.index$alt.no <- alt.no
e.snp.index$all.reads <- e.snp.index$ref.no+e.snp.index$alt.no
e.snp.index$index <- e.snp.index$alt.no/e.snp.index$all.reads

###去掉总reads小于3大于80的位点
e.snp.index <- e.snp.index[e.snp.index$all.reads>3&e.snp.index$all.reads<80,]

得到染色体window的结果，1Mb大小，100kb step步长

max.pos <- c()
for(i in 1:12){  
    max.pos.chr <- max(e.snp.index$POS[e.snp.index$CHROM==i])
    max.pos <- c(max.pos,max.pos.chr)}

wind.str <- c()
wind.end <- c()
wind.chr <- c()
for(i in 1:12){  
    wind.str.perm <- seq(0, max.pos[i]+1000000,100000)      
    wind.end.perm <-  wind.str.perm + 1000000  
    wind.str <- c(wind.str,wind.str.perm)  
    wind.end <- c(wind.end,wind.end.perm)
    wind.chr.perm <- rep(i,length(wind.str.perm))
    wind.chr <- c(wind.chr,wind.chr.perm)
}
wind.bin <-  data.frame(chr=wind.chr,str=wind.str,end=wind.end)
head(wind.bin)

得到SNP frequency ratio

wind.bin$E.index<-  apply(wind.bin,1,function(x){
   wind.bin.index.E <- e.snp.index[e.snp.index$CHROM==x[1]&e.snp.index$POS>=x[2]&e.snp.index$POS<=x[3],]  
   if(dim(wind.bin.index.E)[1] < 10){return(NA)}
   else{return(sum(wind.bin.index.E$alt.no,na.rm=T)/sum(wind.bin.index.E$ref.no,na.rm=T))}})

PS:野生型表型池，代码类似，就不重复了

plot(wind.bin$E.index/wind.bin$W.index~c(1:8146),pch=20,xaxt='n')
##### 主要结果，到此结束。下面的代码是为了画出染色体的边界并进行标注 #######
chr.win.no <- cumsum(table(wind.bin$chr))[1:11]
chr.win.id <- cumsum(table(wind.bin$chr)) - 0.5*table(wind.bin$chr)
for(i in chr.win.no){abline(v=i,lty=2)}
axis(1,at=chr.win.id,label=paste('chr',1:12,sep=''))

原文结果：

注意：
1.细心的同学会发现我的结果虽然与原文类似，但是还是有很大差别的。原因在于原文中定位用到的两个亲本，都不是参考基因组的品种材料，但是野生型和参考基因组类似，我这里为了方便，在统计野生型等位基因频率和突变型等位基因频率时候，用ref reads 代替了wild reads ；用alt reads 代替了mutation reads，仅为演示。
正确的做法应该是，把P1和P2（就是突变型和野生型）分别测序，call出来SNP，用两个亲本间的SNP做后续分析。至于等位基因频率的计算，在下不才，只用了R的代码，R处理起来还是很慢的，大家也可以用perl python等文本处理能力更强大的代码提取出来有效信息，然后再用R做定位分析和画图。当然语言只是工具，只要明白了原理，知道怎么去抽取有效信息，具体怎么出来相信对大家都不是难事。

2.本文用到的samtools bcftools版本相对比较低，更新到最新版本后，bcftools有了质量过滤选项，根据文章后面提供的信息给出如下代码供参考（具体参数含义，请参考官方文档）：

samtools mpileup  -R -d 1000000 -t DP,AN -Q 0 -Bug -f  tomato.SL2.50.fa ejmt.sorted.redup.bam | bcftools  call -mv -O z > ejmt.samtools.raw.vcf
bcftools filter -g 100 -e 'MIN(MQ<50)' ejmt.samtools.raw.vcf | bcftools view -m 2 -M 2 -v snps -o  ejmt.samtools.filtered.vcf

欢迎大家拍砖，提意见。如果大家觉得写的还可以，也请多鼓励，下次聊聊BSA课题设计中需要注意的问题。