接著上一篇留下的問題怜跑，如何按照以1M為窗口沼沈，每次移動(dòng)10Kb計(jì)算均值单绑。我們以KY131是"0/0", 而 "DN422" 是 "1/1"的位點(diǎn)結(jié)果為例傲武，我們來增加這條線
原本我的想法是使用split根據(jù)每個(gè)位點(diǎn)的位置信息進(jìn)行拆分蓉驹，但是這里有一個(gè)問題是城榛，如果只是每個(gè)窗口沒有重疊的話，那么每個(gè)位點(diǎn)應(yīng)該都有一個(gè)唯一的因子戒幔，但這里是以10kb將長度為1Mb的窗口從左往右移動(dòng)吠谢，也就意味著有很多位點(diǎn)時(shí)屬于多個(gè)窗口的土童，也就是無法直接用split-apply-combine的方法進(jìn)行處理。
于是乎，我寫了一個(gè)函數(shù)用于根據(jù)窗口來計(jì)算均值缴淋，它的邏輯就是先創(chuàng)建好窗口资厉，然后根據(jù)窗口的起始位置和終止位置去篩選符合要求的位點(diǎn)，最后進(jìn)行求和罢吃。同時(shí)考慮有些位置可能沒有SNP楚午，那么最后還進(jìn)行了一波過濾。
函數(shù)的輸入是之前snp-index的位置和對(duì)應(yīng)的值尿招，以及確定窗口的大小和步長矾柜，
calcValueByWindow <- function(pos, value,
window_size = 1000000,
step_size = 100000){
# get the max position in the postion
max_pos <- max(pos)
# construct the window
window_start <- seq(0, max_pos + window_size,step_size)
window_end <- window_start + step_size
mean_value <- vector(mode = "numeric", length = length(window_start))
# select the value inside the window
for (j in seq_along(window_start)){
pos_in_window <- which(pos > window_start[j] &
pos < window_end[j])
value_in_window <- value[pos_in_window]
mean_value[j] <- mean(value_in_window)
}
# remove the Not A Number position
nan_pos <- is.nan(mean_value)
mean_value <- mean_value[! nan_pos]
window_pos <- ((window_start + window_end)/ 2)[!nan_pos]
df <- data.frame(pos = window_pos,
value = mean_value)
return(df)
}
得到結(jié)果就可以用lines在之前結(jié)果上加上均值線
par(mfrow = c(3,4))
for (i in paste0("chr", formatC(1:12, width = 2, flag=0)) ){
freq_flt <- freq2[grepl(i,row.names(freq2)), ]
pos <- as.numeric(substring(row.names(freq_flt), 7))
snp_index <- freq_flt[,1] - freq_flt[,2]
# bin
df <- calcValueByWindow(pos = pos, value = snp_index)
plot(x = pos, y =snp_index,
ylim = c(-1,1),
pch = 20, cex = 0.2,
xlab = i,
ylab = expression(paste(Delta, " " ,"SNP index")))
lines(x = df$pos, y = df$value, col = "red")
}
最后再對(duì)文章總結(jié)一下。文章并不是只用了BSA的方法進(jìn)行定位就谜，他們花了幾年的時(shí)間用SSR分子標(biāo)記確定了候選基因可能區(qū)間怪蔑，用BSA的方法在原有基礎(chǔ)上縮小了定位區(qū)間。當(dāng)然即便如此丧荐，候選基因也有上百個(gè)缆瓣，作者通過BLAST的方式，對(duì)這些基因進(jìn)行了注釋虹统。盡管中間還加了一些GO富集分析的內(nèi)容弓坞，說這些基因富集在某個(gè)詞條里，有一個(gè)是DNA metabolic processes(GO:0006259)车荔，但我覺得如果作者用clusterProfiler做富集分析渡冻，它肯定無法得到任何富集結(jié)果。他做富集分析的目的是其實(shí)下面這個(gè)描述忧便，也就是找到和抗凍相關(guān)的基因
LOC_Os06g39740 and LOC_Os06g39750,were annotated as the function of “response to cold (GO: 0009409)”, suggesting their key roles in regulating cold tolerance in rice. "
當(dāng)然他還做了qRT-PCR進(jìn)行了驗(yàn)證族吻，最后推測LOC_Os06g39750應(yīng)該是目標(biāo)基因，這個(gè)基因里還有8個(gè)SNP位點(diǎn)茬腿。
版權(quán)聲明：本博客所有文章除特別聲明外呼奢，均采用 知識(shí)共享署名-非商業(yè)性使用-禁止演繹 4.0 國際許可協(xié)議 (CC BY-NC-ND 4.0) 進(jìn)行許可。
