甲基化数据分析

R脚本-生成密度分布表格

2020-09-09  本文已影响0人  ShawnMagic

在组学分析中，我们常常会看到曼哈顿图，密度分布图，bin图等，之前我以为ggplot2画density plot的数据处理格式和SNP密度一样，后来发现不是那么麻烦。或许是搜索姿势不对，没有找到现成的R包来做这个事情。索性就自己写了个，这个用途应该不局限于密度表格那么简单，不过鉴于现在分析的项目，基本给定一个染色体位置信息，染色体长度信息，自己再设定一个bin的大小，通过这个脚本基本可以实现一键出表格。但是由于用的只是dplyr和Rbase包，通过对数据框的处理而实现的，速度肯定会比较慢。如果在处理大数据比如snp密度等，可能会比较慢。现在你可以用他做：

准备数据

这里以计算陆地棉基因密度分布为例，首先准备好一个下面的表格，需要两列，一列是染色体的名称，另一列是基因的起始位置。

chr start
A01 45486
A01 60648
A01 66475
A01 73986
A01 80100
A01 84484
... ...

此外我们还需要知道每个染色体的长度，这个用TBtools的GXF工具中的获得基因位置和染色体长度的工具就可以得到了，如下：

A01 115949770
A02 105669191
A03 110115461
A04 84614990
A05 108814984
A06 126075219
A07 96696994
A08 125568196
A09 83540628

这样我们就得到了需要的信息。思路如下：
首先我们先按照bin的大小把染色体分成等宽的区块，然后再获得每个区块内基因的个数，那么最重要的两个参数：
binsize，maxVal
想清楚这些后就开始写脚本，如下：

脚本文件

suppressMessages(library(dplyr))
suppressMessages(library(getopt))
options(stringsAsFactors = F)
command=matrix(c(
  'help', 'h', 0, 'logic', 'help information',
  'data', 'd', 1, 'character', 'inputfile: table with headers, the 1st column is label,the second is location',
  'binsize', 'b', 1, 'integer', 'Bin Size: the bin size, if you want set the bin size as 10kb ==>  eg:10000',
  'key', 'k', 1, 'character', 'Key: the lab name you want to extract,if you want calculate chr A01 ==> , eg: "A01"',
  'maxVal', 'm', 1, 'integer', 'Max value : the max value of selected key, such as the length of A01 is 115949770'
),byrow = T, ncol = 5)
args = getopt(command)
## help information
if (!is.null(args$help)) {
  cat(paste(getopt(command, usage = T), "\n"))
  #  q(status=1)
}

## name values
x <- args$data
binsize <- args$binsize
key = args$key
maxVal = args$maxVal
## import data
data = read.table(x,header = T,
                  sep = "\t")
# function
getdensity = function(data,binsize,key,maxVal){
  tmp1 = data.frame(filter(data,data[,1] == key),
                    count = 1)
  bin = seq(binsize,maxVal,by = binsize)
  tmp2 = data.frame(key = key,
                    bin = bin,
                    count = 0)
  for (i in 1:length(bin)) {
    x = (filter(tmp1,tmp1[,2]>bin[i]-binsize & tmp1[,2]<bin[i])[,-1]%>%apply(.,2,sum))[2]
    tmp2[i,3] =  x
  }
  tmp2[i+1,1] = key
  tmp2[i+1,2] = maxVal
  tmp2[i+1,3] = (filter(tmp1,tmp1[,2]>bin[i])[,-1]%>%apply(.,2,sum))[2]
  return(tmp2)
}
## excute
out = getdensity(data = data,
                 binsize = binsize,
                 key = key,
                 maxVal = maxVal)
write.table(out,paste(key,"_",binsize,".Dens.xls",sep = ""),
            row.names = F,
            col.names = F,
            sep = "\t",
            quote = F)

脚本解读

使用方法

Rstudio 中

方法就比较简单了，可以在Rstudio中调整运行，用下面的代码设置相应的参数，比如这里提取的是 A01 的， key 设置成"A01", binsize 设置成10000， maxVal 设置成115949770

suppressMessages(library(dplyr))
suppressMessages(library(getopt))
## import data
data = read.table(x,header = T,
                  sep = "\t")
# function
getdensity = function(data,binsize,key,maxVal){
  tmp1 = data.frame(filter(data,data[,1] == key),
                    count = 1)
  bin = seq(binsize,maxVal,by = binsize)
  tmp2 = data.frame(key = key,
                    bin = bin,
                    count = 0)
  for (i in 1:length(bin)) {
    x = (filter(tmp1,tmp1[,2]>bin[i]-binsize & tmp1[,2]<bin[i])[,-1]%>%apply(.,2,sum))[2]
    tmp2[i,3] =  x
  }
  tmp2[i+1,1] = key
  tmp2[i+1,2] = maxVal
  tmp2[i+1,3] = (filter(tmp1,tmp1[,2]>bin[i])[,-1]%>%apply(.,2,sum))[2]
  return(tmp2)
}
## excute
data = data
binsize = 100000
key = "A01"
maxVal = 115949770
out = getdensity(data = data,
                 binsize = binsize,
                 key = key,
                 maxVal = maxVal)
write.table(out,paste(key,"_",binsize,".Dens.xls",sep = ""),
            row.names = F,
            col.names = F,
            sep = "\t",
            quote = F)

这个后续可以用ggplot2绘图，或者其他处理。

bash下运行

bash下面把第一个代码保存成一个叫densityTable.R的文件，确保你在Rstudio下正确安装了getopt和dplyr两个包后运行

## 如果没有安装getopt包和dplyr包，先打开Rstudio运行：
install.packages("getopt")
install.packages("dplyr")
## 安装成功没有报错后把准备好的第一个文件和这个脚本放到同一个路径下（文件夹中）
## 比如包含染色体编号和基因起始位置的文件叫Gh.gene.pos.xls
Rscript densityTable.R -d Gh.gene.pos.xls -b 100000 -k "A01" -m 115949770
## 运行完你会发现目录中多了一个A01_100000.Dens.xls

结果如下

A01 100000 8
A01 200000 20
A01 300000 11
A01 400000 9
A01 500000 13
A01 600000 8
A01 700000 13
A01 800000 15
A01 900000 11
A01 1000000 19

这样我们就得到了染色体位置，bin位置和基因数目的一个表格，这个表可以用circos画密度分布，热图等等。

批量使用

如果熟悉R的可以直接在R里面用for循环完成，既然讲一键，这里就讲下在bash下面更加简单的while循环出图把
首先看脚本。

while read id
do
        key=$(echo $id | cut -d" " -f 1)
        maxVal=$(echo $id | cut -d" " -f 2)
        echo ${key}
        Rscript densityTable.R -d $2 -b $3 -k ${key} -m ${maxVal}
done <$1
cat *Dens* > Gh.gene_density.txt
rm *Dens*

## Gh.chrlen.xls就是最上面第二个表格，是染色体编号对染色体长度。
## Gh.gene.pos.txt就是第一个表格，是染色体编号对基因起始位置。
sh getgenedensity.sh Gh.chrlen.xls Gh.gene.pos.txt 100000
上一篇 下一篇

猜你喜欢

热点阅读