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')
return(ClusterCut)
}