使用QTLseqr进行BSA-seq分析

QTLseqr是一个用于BSA-seq的R包,它能够读取GATK输出结果进行过滤,最终使用SNP-index和G 统计值的方法进行定位。我们这次尝试使用这个R包来替代我们自己手撸代码来分析一批水稻BSA数据,文章的数据在如何使用BSA方法进行遗传定位(水稻篇)提供了下载链接。

首先是安装和加载R包,由于QTLseqr目前托管在GitHub上,因此需要用devtools才能安装。

devtools::install_github("bmansfeld/QTLseqr")
library("QTLseqr")

在读取数据上,QTLseqr提供了两个函数,importFromGATK和importFromTable,而且这两个函数都没有设置亲本样本名的参数,只有highBulk和lowBulk设置两个混池的样本名。这就很尴尬了,因为我发现直接用importFromGATK读取有4个样本的VCF文件会报错。

不过没关系,我们先将数据预处理成importFromTable能够识别的格式,就可以进行读取了。

我们先用vcfR提取所需要的信息列,包括AD和GT,前者是等位基因深度,后者是基因型

library(vcfR)
vcf <- read.vcfR("04-variant-filter/snps.vcf.gz")

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")

补充下masplit函数的用法,这函数来自于vcfR, 他的作用就是对数据进行分列,注意, 它返回的是矩阵中每个元素拆分后其中一个值,主要参数为

  • myMat 输入矩阵
  • delim = "," 矩阵里元素的分隔符号,
  • count = 0L: 0或1,也就是False或True,统计每个元素有多少个子元素,例如'1,2,3' 会返回3, '1,2' 返回2, 一般用不到
  • record = 1L: 返回分隔之后的第几个元素,
  • sort = 1L: 是否排序
  • decreasing = 1L: 是否降序

之后构建数据框, 过滤亲本杂合的位点,并保存到本地

df <- data.frame(CHROM = chrom,
                 POS = pos,
                 REF = ref,
                 ALT = alt,
                 AD_REF.SRR6327817 = ref_split[,3],
                 AD_ALT.SRR6327817 = alt_split[,3],
                 AD_REF.SRR6327818 = ref_split[,4],
                 AD_ALT.SRR6327818 = alt_split[,4]
                 )

mask <- which(gt[,"SRR6327815"] != "0/1" &  gt[,"SRR6327816"] == "0/1")

df <- df[mask,]
write.table(df, file = "rice.tsv", sep = "\t", row.names = F, quote = F)

之后用importFromTable读取,并过滤掉SNPindex为NaN的行

df <- importFromTable("rice.tsv",
                     highBulk = "SRR6327817",
                     lowBulk = "SRR6327818",
                     chromList = paste0("chr", formatC(1:12, width = 2, flag=0)),
                     sep = "\t")

df <- subset(df, !is.na(SNPindex.LOW) & !is.na(SNPindex.HIGH))

分别运行两种BSA的分析方法,windowSize指的是窗口大小

df <- runGprimeAnalysis(SNPset = df,
                             windowSize = 1e6,
                             outlierFilter = "deltaSNP")

df <- runQTLseqAnalysis(SNPset = df,
                        windowSize = 1e6,
                        popStruc = "RIL",
                        bulkSize = c(20,20))

最后我们可视化结果

plotQTLStats(
  SNPset = df,
  var = "Gprime",
  plotThreshold = TRUE,
  q = 0.01
)

G统计值
plotQTLStats(
  SNPset = df,
  var = "deltaSNP",
  plotIntervals = TRUE)
snp-index-plot

结果和文章比较吻合,在6号染色体有一个非常明显的峰。

QTLseqr这个R包整体的表现还是可以的,就是对输入有一定的要求,如果直接提供VCF文件需要保证你的VCF文件只有两个样本。也就是你需要先单独对亲本做变异检测,然后用这些位点对两个混池变异检测结果进行过滤,不然这两个亲本的简直就体现不出来了。


版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

扫码即刻交流
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 155,704评论 4 358
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,084评论 1 286
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 105,560评论 0 237
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,360评论 0 201
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,686评论 3 285
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,150评论 1 204
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,528评论 2 306
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,256评论 0 193
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,878评论 1 235
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,215评论 2 239
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,757评论 1 255
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,121评论 2 249
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,659评论 3 228
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,932评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,638评论 0 191
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,159评论 2 263
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,112评论 2 258

推荐阅读更多精彩内容