I'm trying to count the GC in regions leading up to polyA clusters. I found the bedtools command nuc, one of the outputs of which is the GC content.

I took the mouse polyA cluster database from here, and the mm10 primary assembly soft masked from here. I imported the polyA clusters bed file to R. I then got a list of all the polyA clusters in the genes I'm interested in. For positive strands on polyA I took the start position, took 100 away from it and made that the start of the region of interest, then made the start of the polyA signal the end. For things on the negative strand I added 100 to the end of the polyA signal. I then exported the results to a bed file format. The code I used to do this is below.

When I try to run:

bedtools nuc -fi Mus_musculus.GRCm38.dna_sm.primary_assembly.fa -bed PolyA.bed

Every signal returns some variant of: Feature (chr10:21,080,315-21,080,415) beyond the length of chr10 size (0 bp). Skipping.

I can however look at several of these regions in the UCSC genome browser without a problem, at least I don't think I had a problem.

Can anyone point me to what's going wrong or to another way I can solve this.

GetRanges <- function(FilteredClusters){
  for (x in 1:nrow(FilteredClusters)){
    if(FilteredClusters$Strand[x] == '+'){
  FilteredClusters$RegionEnd[x] <- FilteredClusters$ClusterStart[x]
  FilteredClusters$RegionStart[x] <- FilteredClusters$RegionEnd[x] - 100
else if(FilteredClusters$Strand[x] == '-')
      FilteredClusters$RegionStart[x] <- FilteredClusters$ClusterEnd[x]
      FilteredClusters$RegionEnd[x] <- FilteredClusters$RegionStart[x] + 100
    ClusterCut <- select(FilteredClusters, 'Chr', 'RegionStart', 'RegionEnd')


There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.

| improve this answer | |
  • $\begingroup$ Well...that is depressingly simple. Thanks for the tip. It's worked just fine. $\endgroup$ – Sethzard May 29 '19 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.