CODEWITHSUNDEEP
r
user3816990
user3816990

Reputation: 247

loop through a file in R and select values based on a range

I have a file called 'genes' with 4300 lines that looks like

Gene_id  chr  start   stop              
GeneA chr1  10  1000                 
GeneB chr1  2300  7000                     
GeneC chr1 10000 13577

and another file called 'bases' (~100,000 lines) and looks like

Chr Bases          
chr1 160           
chr1 157             
chr1 8500           
chr1 2200

I want to produce a file that keeps the bases that fall into the range between start and stop for each gene

so output would look like

Chr Bases             
chr1 160            
chr1 157

I have tried this function but it only gives me back the first entry four times:

methC <- apply(bases,1,function(a){
my_bases <- bases[bases[1]==genes$chr & bases[2]>=genes$start & bases[2]<=genes$stop,]
result <- my_bases[,2]
return(result)
})

>methC
# 160 160 160

so I'm missing base 157 and 160 is duplicated 4 times.

If I use

b <- bases[which(bases[1]==genes$chr & bases[2]>=genes$start & bases[2]<=genes$stop),]
> b
 #  Chr Bases
#chr1   160

I'm still missing 157, but maybe that is due to the order.

However if I try using my real and much larger files with this 'which' function I get back an empty data.frame

> b
Chr        Base       
<0 rows> (or 0-length row.names)

, so that's why I thought a function would be better to handle large datasets.

Upvotes: 1

Views: 789

Answers (1)

Jonny Jones
Jonny Jones

Reputation: 76

I would use the data table library:

library("data.table")

# Read the data from file
genes <- data.table(read.csv("genes.csv"))
bases <- data.table(read.csv("bases.csv"))

# Define the columns that are used to join the two tables
setkey(genes, chr)
setkey(bases, Chr)

# The first line joins the two tables, and the second line
# filters the result to keep only those that match and defines
# the columns that are to be shown in the output. Note the need
# to take care with capitalisation and the reserved word 'stop'
genes[bases, allow.cartesian = TRUE][
  Bases >= start & Bases <= `stop`, list(Chr = chr, Bases)]

Which produces this result:

    Chr Bases
1: chr1   160
2: chr1   157

Upvotes: 2

Related Questions