BSA分析拟南芥F2代分离群体混池测序

1. 实验背景

为了研究拟南芥对高温响应的基因，我们对拟南芥的野生型Col进行了EMS诱变，通过对诱变后的种子多代的高温筛选，我们发现了一个对高温敏感的突变体，该突变体的下胚轴的长度在高温下要比野生型显著的短。之后，将此突变体和野生型Col进行杂交，F1表现长下胚轴，F1自交，F2出现了明显的性状分离，即表现长下胚轴和短下胚轴两种类型（长：短～3:1），遗传分析表明该突变是一个隐形突变，有单基因控制。

2. 实验设计及测序

对F2群体中的长，短下胚轴的两种类型的材料分别取30株，然后混合提取DNA，建立两个DNA池，long-pool, short-pool。之后选取亲本Col，及突变体进行建库测序。一共四个样品，采用ILUMINA双端测序。每个材料测序40～50X。公司返回的数据，每个样品大约是7Gb.根据拟南芥基因组的大小125Mb，本次测序每个样品的深度大约是56X。返回的原始数据如下：

mkdir BSA_project
cd BSA_project
mkdir Rawdata
#move your raw data here
cd Rawdata/
ls
Cf-long_R1.fq.gz  Cf-short_R1.fq.gz  Col_R1.fq.gz  mutant_R1.fq.gz  md5.txt
Cf-long_R2.fq.gz  Cf-short_R2.fq.gz  Col_R2.fq.gz  mutant_R2.fq.gz

3. 数据分析。

（1）创建序列回帖的参考基因组index， GATK call SNP 的index。根据参考基因组fastq名称运行一下脚本

cd BSA_project
mkdir ref #参考基因组文件，INDEX，GATK的dict等
cd ref
ls
Athaliana_447_TAIR10.fa
mkdir script
# put scripts here

#!/bin/bash
# building sequence alginment dictionary, samtools faidx and gatk creatSequenceDictionary
#Usage: sh gatk_step1.sh /path/your_genome.fasta
bwa=/home/zhanghuairen/bin/bwa                   # set where to find software
gatk=/home/zhanghuairen/software/gatk-4.1.7.0/gatk
samtools=/home/biosoftware/bin/samtools#bwa index
reference=$1
time $bwa index "$reference" && echo "** bwa index done! ** "
#samtools index
time  $samtools faidx $reference && echo "** samtools faidx done! ** "#注意：使用GATK之前，需要先建立参考基因组索引文件.dict和.fai
#.dict中包含了基因组中contigs的名字，也就是一个字典；
#.fai也就是fasta index file，索引文件，可以快速找出参考基因组的碱基，由samtools faidx构建
#构建.dict文件（原来要使用picard的CreateSequenceDictionary模块，但是现在gatk整合了此模块，可以直接使用）
# gatk createSequenceDictionary
time $gatk --java-options "-Xmx100G -Djava.io.tmpdir=./tmp" CreateSequenceDictionary \-R "$reference" \-O "$reference.dict" \&& echo "** gatk createSequenceDictionary done! **"

在上面的ref文件夹中运行该脚本,会生成bwa比对的参考基因组文件的INDEX。以及GATK所需要的dict.这个时候要把GATK的dict 该一个名称，比如：mv Athaliana_447_TAIR10.fa.dict Athaliana_447_TAIR10.dict。不然下边GATKcall SNP 会报错

运行完之后的ref 包含如下：

 Athaliana_447_TAIR10.fa.amb      Athaliana_447_TAIR10.fa.pacAthaliana_447_TAIR10.fa.ann      Athaliana_447_TAIR10.fa.sa
Athaliana_447_TAIR10.dict          Athaliana_447_TAIR10.fa.bwt
Athaliana_447_TAIR10.fa            Athaliana_447_TAIR10.fa.fai

（2）对每个原始数据进行质控，去除接头，低质量的reads; bwa进行比对；samtools 对bam文件进行排序、索引。最后使用GATK call SNP 生成GVCF文件。这个脚本如下：

#!/bin/bash
# this is a whole gonome gatk snp calling full process
#Uage: sh gatk_step2.sh sample_1.fastq sample_2.fastq out_dictionary
fastp=/usr/local/bin/fastp
bwa=/home/zhanghuairen/bin/bwa
gatk=/home/zhanghuairen/software/gatk-4.1.7.0/gatk
samtools=/home/biosoftware/bin/samtools
reference=/home/zhanghuairen/lujingyun/ref/Athaliana_447_TAIR10.fa
NTHREADS=30
fq1=$1
fq2=$2
sample=${fq1%%_*}
outdir=$3
outdir=${outdir}/${sample}if [ ! -d "$outdir/cleanfq" ]
then mkdir -p "$outdir/cleanfq"
fiif [ ! -d "$outdir/bwa" ]
then mkdir -p "$outdir/bwa"
fiif [ ! -d "$outdir/gatk" ]
then mkdir -p "$outdir/gatk"
fitime $fastp -i "$fq1" \-I "$fq2" \-o "$outdir/cleanfq/${sample}_1.fastq.gz" \-O "$outdir/cleanfq/${sample}_2.fastq.gz" \-h "$outdir/cleanfq/${sample}.html" \-j "$outdir/cleanfq/${sample}.json" \&& echo '** fq QC done'
time $bwa mem -M -R "@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA" \-t $NTHREADS \$reference \"$outdir/cleanfq/${sample}_1.fastq.gz"  "$outdir/cleanfq/${sample}_2.fastq.gz" \| $samtools view -bS - > "$outdir/bwa/${sample}.bam"\&& echo '** BWA MEM done ** '
time $samtools sort \-@ $NTHREADS \-O bam \-o "$outdir/bwa/${sample}.sorted.bam" \"$outdir/bwa/${sample}.bam"\&& echo "** sorted raw bamfiles done"
time $samtools index "$outdir/bwa/${sample}.sorted.bam" && echo "** ${sample}.sorted.bam index done "# 去除PCR重复
time $gatk MarkDuplicates \-I "$outdir/bwa/${sample}.sorted.bam" \-M "$outdir/bwa/${sample}.markup_metrics.txt" \-O "$outdir/bwa/${sample}.sorted.markup.bam" \1> "$outdir/bwa/${sample}.log.mark" 2>&1 \&& echo "** ${sample}.sorted.bam MarkDuplicates done **"# samtools index 去除PCR标记的bam 文件
time $samtools index "$outdir/bwa/${sample}.sorted.markup.bam" \&& echo " ** ${sample}.sorted.markup.bam index done **"#gatk开始：必选 -I -O -R，代表输入、输出、参考
#接下来可以按照字母顺序依次写出来，这样比较清晰
#-bamout：将一整套经过gatk程序重新组装的单倍体基因型（haplotypes）输出到文件
#-stand-call-conf :低于这个数字的变异位点被忽略，可以设成标准30（默认是10）
time $gatk --java-options "-Xmx100G -Djava.io.tmpdir=./tmp" HaplotypeCaller \-R $reference \-I "$outdir/bwa/${sample}.sorted.markup.bam" \-O "$outdir/gatk/${sample}.HC.gvcf.gz" \--tmp-dir "$outdir/gatk"\--emit-ref-confidence GVCF \-stand-call-conf 10 \&& echo "** GVCF ${sample}.HC.gvcf.gz done ***"

对每一个样品运行上面的脚本，然后会在输出的文件夹中创建四个样品的文件名，每个文件名之内包含三个文件夹，分别是bwa, cleanfq, gatk.

├── Cf-long
│   ├── bwa
│   │   ├── Cf-long.bam
│   │   ├── Cf-long.log.mark
│   │   ├── Cf-long.markup_metrics.txt
│   │   ├── Cf-long.sorted.bam
│   │   ├── Cf-long.sorted.bam.bai
│   │   ├── Cf-long.sorted.markup.bam
│   │   └── Cf-long.sorted.markup.bam.bai
│   ├── cleanfq
│   │   ├── Cf-long_1.fastq.gz
│   │   ├── Cf-long_2.fastq.gz
│   │   ├── Cf-long.html
│   │   └── Cf-long.json
│   └── gatk
│       ├── Cf-long.HC.gvcf.gz
│       └── Cf-long.HC.gvcf.gz.tbi
├── Cf-short
│   ├── bwa
│   │   ├── Cf-short.bam
│   │   ├── Cf-short.log.mark
│   │   ├── Cf-short.markup_metrics.txt
│   │   ├── Cf-short.sorted.bam
│   │   ├── Cf-short.sorted.bam.bai
│   │   ├── Cf-short.sorted.markup.bam
│   │   └── Cf-short.sorted.markup.bam.bai
│   ├── cleanfq
│   │   ├── Cf-short_1.fastq.gz
│   │   ├── Cf-short_2.fastq.gz
│   │   ├── Cf-short.html
│   │   └── Cf-short.json
│   └── gatk
│       ├── Cf-short.HC.gvcf.gz
│       └── Cf-short.HC.gvcf.gz.tbi
├── Col
│   ├── bwa
│   │   ├── Col.bam
│   │   ├── Col.log.mark
│   │   ├── Col.markup_metrics.txt
│   │   ├── Col.sorted.bam
│   │   ├── Col.sorted.bam.bai
│   │   ├── Col.sorted.markup.bam
│   │   └── Col.sorted.markup.bam.bai
│   ├── cleanfq
│   │   ├── Col_1.fastq.gz
│   │   ├── Col_2.fastq.gz
│   │   ├── Col.html
│   │   └── Col.json
│   └── gatk
│       ├── Col.HC.gvcf.gz
│       └── Col.HC.gvcf.gz.tbi
├── mutant
│   ├── bwa
│   │   ├── mutant.bam
│   │   ├── mutant.log.mark
│   │   ├── mutant.markup_metrics.txt
│   │   ├── mutant.sorted.bam
│   │   ├── mutant.sorted.bam.bai
│   │   ├── mutant.sorted.markup.bam
│   │   └── mutant.sorted.markup.bam.bai
│   ├── cleanfq
│   │   ├── mutant_1.fastq.gz
│   │   ├── mutant_2.fastq.gz
│   │   ├── mutant.html
│   │   └── mutant.json
│   └── gatk
│       ├── mutant.HC.gvcf.gz
│       └── mutant.HC.gvcf.gz.tbi

(3)将上面生成的GVCF文件使用GATK进行整合，得到一个完整的GVCF文件,并使用GenotypeGVCFs 命令整合。运行下边的命令：

mkdir gcvf
# move the single gvcf file into the directory
gatk=/home/zhanghuairen/software/gatk-4.1.7.0/gatk
reference=/home/zhanghuairen/lujingyun/ref/Athaliana_447_TAIR10.fa
cd gcvf
time $gatk CombineGVCFs \-R $reference \-V Cf-long.HC.gvcf.gz \-V Cf-short.HC.gvcf.gz \-V Col.HC.gvcf.gz \-V mutant.HC.gvcf.gz \-O BSA.gvcf.gz
time $gatk GenotypeGVCFs \-R  $reference \-V "BSA.gvcf.gz" \-O "BSA_no_filter.HC.vcf"

(4)至此，我们已经形成了一个初步的vcf 文件，之后就是对这个文件进行过滤，这里我们仅仅对这些SNP和INDEL进行硬过滤，也不适用其他的矫正方法。

#!/bin/bash
gatk=/home/zhanghuairen/software/gatk-4.1.7.0/gatk
reference=/home/zhanghuairen/lujingyun/ref/Athaliana_447_TAIR10.fa
time $gatk SelectVariants -R $reference \-V  BSA_no_filter.HC.vcf \-select-type SNP \-O BSA_combine_SNPs.vcf
#select INDELS
time $gatk SelectVariants -R $reference \-V  BSA_no_filter.HC.vcf \-select-type INDEL  \-O BSA_combine_INDELs.vcf
#To filter SNPs: 其他的一些硬过滤的指标不知道为什么在此报错，我查看了一下VCF文件，可能并不是每个标记都有哪些过滤的指标。
time $gatk VariantFiltration  \-V BSA_combine_SNPs.vcf \-filter "MQ < 30.0" \--filter-name "MQ_filter_SNP" \-filter  "QD < 2.0" \--filter-name "QD_filter_SNP" \-O BSA_SNPs_filter.vcfgrep -E '^#|PASS' BSA_SNPs_filter.vcf > BSA_SNPs_filterPASSED.vcf#filter INDELS
time $gatk VariantFiltration \-V BSA_combine_INDELs.vcf \-filter  "MQ < 30.0" \--filter-name "MQ_filter_INDEL" \-filter "SOR > 10.0" \--filter-name "SOR_filter_INDEL" \-filter "QD < 2.0" \--filter-name "QD_filter_INDEL" \-filter "FS > 200.000" \--filter-name "FS_filter_INDEL" \-O BSA_INDELs_filter.vcf
time $gatk MergeVcfs \。#整合SNP和INDEL-I BSA_INDELs_filterPASSED.vcf \-I BSA_SNPs_filterPASSED.vcf \-O BSA_filter_merged_SNP_INDEL.vcfgrep -E '^#|PASS' BSA_INDELs_filter.vcf > BSA_INDELs_filterPASSED.vcf
#This is the number of sites before filtering:
grep -vc "^#" BSA_combine_SNPs.vcf
grep -vc "^#" BSA_combine_INDELs.vcf
#This is the number of sites retained after filtering:
grep -vc "^#" BSA_SNPs_filter.vcf
grep -vc "^#" BSA_INDELs_filter.vcf#最后我发现上面的过滤完之后，和原来的差别不是很大。具体原因还在学习中。

（5）拿到了过滤后的SNP和INDEL的文件，就可以做BSA的分析了，这里才用R的QTLseqr 包。

setwd("~/BSA_project/R_analysis")
#devtools::install_github("bmansfeld/QTLseqr")
library("QTLseqr")
library(vcfR)
vcf <- read.vcfR("BSA_filter_merged_SNP_INDEL.vcf")
#Cf-long    Cf-short        Col     mutant #前两个分别是两个池，后边的是野生型和突变体
chrom <- getCHROM(vcf)
pos <- getPOS(vcf)
ref <- getREF(vcf)
alt <- getALT(vcf)ad <- extract.gt(vcf, "AD")
ref_split <- masplit(ad, record = 1, sort = 0)
alt_split <- masplit(ad, record = 2, sort = 0)
gt <- extract.gt(vcf, "GT")
head(gt)
df <- data.frame(CHROM = chrom,POS = pos,REF = ref,ALT = alt,AD_REF.long= ref_split[,1],AD_ALT.long= alt_split[,1],AD_REF.short = ref_split[,2],AD_ALT.short = alt_split[,2]
)
mask <- which(gt[,"Col"] != "0/1" &  gt[,"mutant"] != "0/1") # 选择亲本纯合的标记
df_mask <- df[mask,]
write.table(df_mask, file = "BSA_right2.tsv", sep = "\t", row.names = F, quote = F)
df_cal <- importFromTable("BSA_right2.tsv",highBulk = "long",lowBulk = "short",sep = "\t")
df_clean <- subset(df_cal, !is.na(SNPindex.LOW) & !is.na(SNPindex.HIGH))
df_result <- runGprimeAnalysis(SNPset = df_clean,windowSize = 1e6,outlierFilter = "Gprime")
plotQTLStats(SNPset = df_result,var = "Gprime",plotThreshold = TRUE,q = 0.01
)

最后，现在定位结果，定位结果的具体数据可以从df_results 中获得，该文件是一个数据框，详细的记录了每一个标记，不同等位基因的等位基因频率，即SNP index。也包含了G值，P值，以及矫正P值。我们可以根据这些结果清楚的知道突变基因的大概位置。由于群体是突变体和野生型杂交的，有用的标记不是很多，定位区间也很大，但是可以给予个明确的参考。

4致谢

至此，一个比较完整的BSA分析流程就构建了，当然这个流程也并不非常完整，尤其是对SNP和INDEL进行过滤的那一步，我做的还是比较粗糙的。这也是由于没有仔细去学习GATK的软件，以及论坛资料。在本次分析中参考了很多前人写的关于BSA分析的方法，步骤，再次表示衷心感谢。我将以下参考资源列了出来，没有列的只是我忘了，看完网页没有收藏，请多包涵。

https://evodify.com/gatk-in-non-model-organism

/https://www.jianshu.com/p/1b860f137fae

https://zhuanlan.zhihu.com/p/36745259

https://cloud.tencent.com/developer/article/1445600

https://www.jianshu.com/nb/36118798

https://www.jianshu.com/p/2e21e4603457

http://xuzhougeng.top/archives/QTL-analysis-using-BSA-seq-in-rice

https://github.com/caulilin

BSA分析拟南芥F2代分离群体混池测序相关推荐

Mutmap定位拟南芥的基因
mutmap 定位基因也是基于BSA的方法,之前一篇博客 BSA分析拟南芥F2代分离群体混池测序是做的突变体也野生型杂交后代中F2选择极端个体,加亲本进行测序分析的一篇流程.上篇中的数据也可以用mu ...
Nature: 拟南芥微生物组功能研究1培养组学—高通量细菌分离培养鉴定
点击上方蓝色「宏基因组」关注我们! 背景介绍在近年来,随着宏基因组技术的发展,对原核生物的纯培养逐渐被忽视,但纯培养对于阐明这些微生物的功能是必需的. 本文是2015末发表在Nature的Artic ...
Nature：拟南芥微生物组功能研究1培养组学—高通量细菌分离培养鉴定
背景介绍 Bai, Y., et al. (2015). "Functional overlap of the Arabidopsis leaf and root microbiota.&q ...
ISME | 根内生真菌与来自拟南芥和大麦微生物群落协同有益作用
题目:The fungal root endophyte Serendipita vermifera displays inter-kingdom synergistic beneficial eff ...
中国科学：拟南芥二半萜类化合物调控根系微生物组
2019年5月10日关于王国栋组和白洋组二半萜化合物对微生物组调控的文章发表于<中国科学>,新闻稿如下: 中国科学:中科院遗传发育所揭示拟南芥二半萜对根系微生物组的调控机制同一天由白洋组 ...
SCLS：拟南芥二半萜类化合物调控根系微生物组
2019年5月10日关于王国栋组和白洋组二半萜化合物对微生物组调控的文章发表于<中国科学>,新闻稿如下: 中国科学:中科院遗传发育所揭示拟南芥二半萜对根系微生物组的调控机制同一天由白洋组 ...
SCLS：中科院遗传发育在拟南芥二半萜类化合物调控根系微生物组取得突破进展
文章目录摘要结果拟南芥中进化出新的二半萜基因簇图1. 十字花科中广泛共线性的TPS-GFPPS-P450基因簇单个氨酸取代决定二半萜的特异性图2. 十字花科植物二半萜特异的底物拟南芥二半 ...
【Front Plant Sci】BrMYB2调控大白菜和拟南芥发育过程中花青素生物的合成机制
文章信息题目:Dynamic Changes of the Anthocyanin Biosynthesis Mechanism During the Development of Heading ...
Nature：拟南芥微生物组功能研究2细菌基因组测序和分析
本网对Markdown排版支持较差,请跳转植物微生物组公众号阅读背景介绍 Bai, Y., et al. (2015). "Functional overlap of the Arabid ...

BSA分析拟南芥F2代分离群体混池测序

BSA分析拟南芥F2代分离群体混池测序相关推荐

最新文章

热门文章