Tips:ChIPseeker注释的peak文件不要随意命名

2024-10-25  本文已影响0人  生信云笔记

  ChIPseeker用来做peak注释体验感确实很好,所以也经常用,为此还特意包装了一个性化脚本方便命令行直接使用,后台也有不少人获取脚本,省时省力能解决问题的工具当然是不错的选择。

  不过,做为ChIPseeker的老用户,个人觉得在使用过程中有一点需要留意:ChIPseeker对注释文件名是敏感的。免得会丢失第一个peak,虽然少一个应该不会有什么影响,但严谨一点总归是好的。简单来说,就是默认情况下,ChIPseeker会根据文件名的后缀来判断文件是否需要有header。所以,正常情况下,只要不乱给文件添加header和随意命名是不会有问题的,如果喜欢瞎倒腾(比如本人),那就得注意了,读取文件时就要多加一个参数了。

  下面还是用包自带的数据来说明一下,懒得准备数据了,能说明问题就好:

library(ChIPseeker)

peakfile <- getSampleFiles()[[1]]
peakfile
[1] "/home/software/rlib/r-4.0_lib/ChIPseeker/extdata/GEO_sample_data/GSM1174480_ARmo_0M_peaks.bed.gz"
peak <- readPeakFile(peakfile)
peak
GRanges object with 812 ranges and 0 metadata columns:
        seqnames              ranges strand
           <Rle>           <IRanges>  <Rle>
    [1]     chrX   61728297-61728780      *
    [2]    chr10   39105185-39105362      *
    [3]     chrY   13137266-13137499      *
    [4]    chr11 114049918-114050234      *
    [5]     chrY   13107715-13107867      *
    ...      ...                 ...    ...
  [808]     chrX   49239222-49239305      *
  [809]     chrX   54945698-54945789      *
  [810]     chrX   61817143-61817176      *
  [811]     chrX 147048421-147048507      *
  [812]     chrY       887860-887931      *
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

  把上面的示例文件名随意重命名,看看上面的读取方式有什么不同:

file.copy(peakfile, 'GSM1174480_ARmo_0M_peaks.bed3.gz')
peakfile1 <- 'GSM1174480_ARmo_0M_peaks.bed3.gz'
peak1 <- readPeakFile(peakfile1)
peak1
GRanges object with 811 ranges and 0 metadata columns:
        seqnames              ranges strand
           <Rle>           <IRanges>  <Rle>
    [1]    chr10   39105185-39105362      *
    [2]     chrY   13137266-13137499      *
    [3]    chr11 114049918-114050234      *
    [4]     chrY   13107715-13107867      *
    [5]     chr1 142655900-142656170      *
    ...      ...                 ...    ...
  [807]     chrX   49239222-49239305      *
  [808]     chrX   54945698-54945789      *
  [809]     chrX   61817143-61817176      *
  [810]     chrX 147048421-147048507      *
  [811]     chrY       887860-887931      *
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

  同样的文件,文件名不同,用同样的方式读取结果有了一行的差别。上面的第一行在下面丢失了,这种情况应该加一个参数,正确的读取方式如下:

peak1 <- readPeakFile(peakfile1, header=F)
peak1
GRanges object with 812 ranges and 0 metadata columns:
        seqnames              ranges strand
           <Rle>           <IRanges>  <Rle>
    [1]     chrX   61728297-61728780      *
    [2]    chr10   39105185-39105362      *
    [3]     chrY   13137266-13137499      *
    [4]    chr11 114049918-114050234      *
    [5]     chrY   13107715-13107867      *
    ...      ...                 ...    ...
  [808]     chrX   49239222-49239305      *
  [809]     chrX   54945698-54945789      *
  [810]     chrX   61817143-61817176      *
  [811]     chrX 147048421-147048507      *
  [812]     chrY       887860-887931      *
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

  ChIPseeker默认自动判断文件名是否符合无header条件,上面的第二种情况下,默认文件第一行为header所以丢失了。

  反过来也是一样,如果命名符合有header条件,实际内容里却没有,那么也需要添加header=F参数,否则第一行也会做为header而丢失。

  readPeakFile读取文件时,文件名以.bed.bed.gzPeak.gz.bedGraph.gz.narrowPeak.broadPeak.gappedPeak这些后缀的内容默认无header,其他情况默认有header

  通常情况下,readPeakFile除了读取的文件,无需额外提供参数,查看函数的帮忙信息:

Usage:
     readPeakFile(peakfile, as = "GRanges", ...)

Arguments:
peakfile: peak file
      as: output format, one of GRanges or data.frame
     ...: additional parameter

  可知,除了前面确定的两个参数外,还可以提供其他额外的参数,除了这里说的header,其余额外的参数还可以是哪些,可以查看read.delim的帮助信息。

上一篇 下一篇

猜你喜欢

热点阅读