Hi,

I'm new here and wanted to know if it is possible to subset a BSgenome and keep only the standard chromosomes part (1 to 22, X, Y and MT for BSgenome.Hsapiens.NCBI.GRCh38). The goal is also to keep a BSgenome as a result. I use vmatchPattern and vmatchPDict to do some statistics and try to test visual results with ggplot2 library. It would be easier to get only the results I'm interested in. I tried to use seqlevels but you can't because the supplied 'seqinfo' must have the same length as the current 'seqinfo' when replacing the 'seqinfo' If you tell me RTFM I would do it but if it's not easy I appreciate to learn the trick. Thanks a lot.

Alain ```

There is an 'exclude' argument go both vmatchPattern and vmatchPDict that you can use to exclude the patches and alternative haplotypes. Have you tried that? An alternative would be to get the FASTA for GRCh38 and generate your own BSgenome package with just the standard chromosomes.

> library(BSgenome.Hsapiens.UCSC.hg38)
 > data(HNF4alpha)
 > HNF4alpha
 DNAStringSet object of length 71:
      width seq
  [1]    13 AGTTCAAGGATCA
  [2]    13 GGGGTCAAGGGTT
  [3]    13 AGGGTAAAGGTTG
  [4]    13 GTCACAAAAGTCC
  [5]    13 GGTCCAAAGGGCG
  ...   ... ...
 [67]    13 CTTGGAACCGGGG
 [68]    13 AGGTCAGGGTCCC
 [69]    13 TGTCCAAAGTCCA
 [70]    13 TGATCAGACAAAG
 [71]    13 AAACCAAAGTTCA
 > pattern <- consensusString(HNF4alpha)
 > pattern
 [1] "RGKBCAAAGKYCA"
 > vmatchPattern(pattern, Hsapiens, fixed = "subject")
 GRanges object with 7814 ranges and 0 metadata columns:
                      seqnames          ranges strand
                         <Rle>       <IRanges>  <Rle>
      [1]                 chr1     30272-30284      +
      [2]                 chr1 1307002-1307014      +
      [3]                 chr1 1552387-1552399      +
      [4]                 chr1 2022905-2022917      +
      [5]                 chr1 2393725-2393737      +
      ...                  ...             ...    ...
   [7810] chr22_KQ458388v1_alt   104632-104644      -
   [7811] chr22_KQ759761v1_alt       8663-8675      -
   [7812] chr22_KQ759761v1_alt     22332-22344      -
   [7813] chr22_KQ759761v1_alt     75045-75057      -
   [7814]  chrX_KV766199v1_alt   145613-145625      -
   -------
   seqinfo: 640 sequences (1 circular) from hg38 genome

> vmatchPattern(pattern, Hsapiens, fixed = "subject", exclude = seqnames(Hsapiens)[-(1:24)])
 GRanges object with 7316 ranges and 0 metadata columns:
          seqnames            ranges strand
             <Rle>         <IRanges>  <Rle>
      [1]     chr1       30272-30284      +
      [2]     chr1   1307002-1307014      +
      [3]     chr1   1552387-1552399      +
      [4]     chr1   2022905-2022917      +
      [5]     chr1   2393725-2393737      +
      ...      ...               ...    ...
   [7312]     chrY 22554712-22554724      -
   [7313]     chrY 23129977-23129989      -
   [7314]     chrY 24763691-24763703      -
   [7315]     chrY 25094790-25094802      -
   [7316]     chrY 26029106-26029118      -
   -------
   seqinfo: 640 sequences (1 circular) from hg38 genome

I think I have a solution reading the manual. By using bsapply function on a BSgenome object with BSParams and the exclude parameter.

Example :

library("BSgenome.Hsapiens.NCBI.GRCh38")

param <- new("BSParams", X = Hsapiens, FUN = letterFrequency, exclude = "_")

GC.freq <- bsapply(param, "GC", as.prob = T)

Do someone agree with me ? Thanks for telling.