Comparing row values in pandas dataframe

Question

I have data in a pandas dataframe where two columns contain numerical sequences (start and stop). I want to identify which rows have stop values which overlap with the next rows' start values. Then I need to concatenate them into a single row so that I only have single none-overlapping numerical sequences represented by my start and stop values in each row.

I have loaded my data into a pandas dataframe:

  chr     start       stop        geneID
0 chr13   32889584    32889814    BRCA2
1 chr13   32890536    32890737    BRCA2
2 chr13   32893194    32893307    BRCA2
3 chr13   32893282    32893400    BRCA2
4 chr13   32893363    32893466    BRCA2
5 chr13   32899127    32899242    BRCA2

I want to compare the rows in the dataframe. Check whether the stop value for each row is less than the start value for the following row and then create a row in a new dataframe with the correct start and stop values. Ideally when there are several rows which all overlap this would be concatenated all in one go, however I suspect I will have to iterate over my output until this doesn't happen any more.

My code so far can identify whether there is an overlap (adapted from this post):

import pandas as pd
import numpy as np

columns = ['chr','start','stop','geneID']
bed = pd.read_table('bedfile.txt',sep='\s',names=['chr','start','stop','geneID'],engine='python')

def bed_prepare(inp_bed):
    inp_bed['next_start'] = inp_bed['start'].shift(periods=-1)
    inp_bed['distance_to_next'] = inp_bed['next_start'] - inp_bed['stop']
    inp_bed['next_region_overlap'] = inp_bed['next_start'] < inp_bed['stop']
    intermediate_bed = inp_bed
    return intermediate_bed

And this gives me output like this:

print bed_prepare(bed)

       chr     start      stop geneID  next_start  distance_to_next  next_region_overlap
0    chr13  32889584  32889814  BRCA2    32890536               722  False
1    chr13  32890536  32890737  BRCA2    32893194              2457  False
2    chr13  32893194  32893307  BRCA2    32893282               -25  True
3    chr13  32893282  32893400  BRCA2    32893363               -37  True
4    chr13  32893363  32893466  BRCA2    32899127              5661  False

I want to put this intermediate dataframe into the following function in order get the desired output (shown below):

new_bed = pd.DataFrame(data=np.zeros((0,len(columns))),columns=columns)

def bed_collapse(intermediate_bed, new_bed,columns=columns):
    for row in bed.itertuples():
    output = {}
        if row[7] == False:
            # If row doesn't overlap next row, insert into new dataframe unchanged.                                   
            output_row = list(row[1:5])
        if row[7] == True:
            # For overlapping rows take the chromosome and start coordinate                                           
            output_row = list(row[1:3])
            # Iterate to next row                                                                                     
            bed.itertuples().next()
            # Append stop coordinate and geneID                                                                       
            output_row.append(row[3])
        output_row.append(row[4])
        #print output_row                                                                                             
        for k, v in zip(columns,output_row): otpt[k] = v
        #print output                                                                                                 
        new_bed = new_bed.append(otpt,ignore_index=True)
    output_bed = new_bed
    return output_bed


int_bed = bed_prepare(bed)
print bed_collapse(int_bed,new_bed)

Desired output:

  chr     start       stop        geneID
0 chr13   32889584    32889814    BRCA2
1 chr13   32890536    32890737    BRCA2
2 chr13   32893194    32893466    BRCA2
5 chr13   32899127    32899242    BRCA2

However, when I run the function I get my original dataframe back unchanged. I know that the problem is when I try to call bed.itertuples().next(), as this is clearly not quite the right syntax/location for the call. But I don't know the correct way to rectify this.

Some pointers would be great.

SB :)

Update

This is a BED file where each row refers to an amplicon (genomic region) with start and stop coordinates. Some of the amplicons overlap; ie the start coordinate is before the stop coordinate on the previous row. Therefore I need to identify which rows overlap and concatenate the correct starts and stops so that each row represents and entirely unique amplicon which doesn't overlap any other row.

Answer 1

pyranges will allow you to do this super-quickly in one line of code:

import pyranges as pr

c = """Chromosome     Start       End        geneID
chr13   32889584    32889814    BRCA2
chr13   32890536    32890737    BRCA2
chr13   32893194    32893307    BRCA2
chr13   32893282    32893400    BRCA2
chr13   32893363    32893466    BRCA2
chr13   32899127    32899242    BRCA2"""

gr = pr.from_string(c)
# +--------------+-----------+-----------+------------+
# | Chromosome   |     Start |       End | geneID     |
# | (category)   |   (int32) |   (int32) | (object)   |
# |--------------+-----------+-----------+------------|
# | chr13        |  32889584 |  32889814 | BRCA2      |
# | chr13        |  32890536 |  32890737 | BRCA2      |
# | chr13        |  32893194 |  32893307 | BRCA2      |
# | chr13        |  32893282 |  32893400 | BRCA2      |
# | chr13        |  32893363 |  32893466 | BRCA2      |
# | chr13        |  32899127 |  32899242 | BRCA2      |
# +--------------+-----------+-----------+------------+
# Unstranded PyRanges object has 6 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

m = gr.merge(by="geneID")
# +--------------+-----------+-----------+------------+
# | Chromosome   |     Start |       End | geneID     |
# | (category)   |   (int32) |   (int32) | (object)   |
# |--------------+-----------+-----------+------------|
# | chr13        |  32889584 |  32889814 | BRCA2      |
# | chr13        |  32890536 |  32890737 | BRCA2      |
# | chr13        |  32893194 |  32893466 | BRCA2      |
# | chr13        |  32899127 |  32899242 | BRCA2      |
# +--------------+-----------+-----------+------------+
# Unstranded PyRanges object has 4 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Note that by="geneID" makes it so intervals are only merged if they overlap and have the same value for geneID. See also the method cluster if you want to merge the interval meta-data with a custom function.

Answer 2

I am not sure I understand why you are doing what you are doing, but you can get your desired output by simply using indexing. e.g.

# assume your data is stored in <df>
# call the temporary dataframe <tmp>
tmp = df[ ['chr','start','stop','geneID'] ][(df.stop - df.start.shift(-1))>0]

Is that what you are trying to do, ultimately?

UPDATE Ok, I see what you are doing. Bear in mind that I have never worked with any genome data, so I have no idea how many rows are in your columns so simple "looping" may be quite slow (if you have a few billion rows this could take a while), but it is the only solution that comes to mind. Here is the first thing to come to mind (NOTE: this is not a finished product since you need to determine how to handle the NaN's that are introduced and how to handle the loop termination).

import pandas as pd

df = pd.DataFrame(index = [0,1,2,3,4,5],columns=['chr','start','stop','geneID'])

df['chr']    = np.array( ['chr13']*6 )
df['start']  = np.array( [32889584,32890536,32893194,32893282,32893363,32899127] )
df['stop']   = np.array( [32889814,32890737,32893307,32893400,32893466,32899242] )
df['geneID'] = np.array( ['BRCA2']*6 )

# calculate difference between start/stop times for adjacent rows
# this will effectively "look into the future" to see if the upcoming row has 
# a start time that is greater than the current stop time
df['tdiff'] = (df.start - df.stop.shift(1)).shift(-1)

# create new dataframe
df_cut = df.copy()*0

r = 0
while r < df.shape[0]:
    if df.tdiff[r] > 0:
        df_cut.iloc[r] = df.iloc[r]
        r+=1

    elif df.tdiff.iloc[r] < 0: # have to determine how you will handle the NaN's later
        df_cut.chr.iloc[r] = df.chr.iloc[r]
        df_cut.start.iloc[r] = df.start.iloc[r]
        df_cut.geneID.iloc[r] = df.geneID.iloc[r]

        # get the next-valid row and put "stop" value into <df_cut>
        df_cut.stop.iloc[r] = df.ix[r:][df.tdiff>0].stop.iloc[0]

        # determine new index location for <r>
        r = df.ix[r:][df.tdiff>0].index[0] + 1

# eliminate empty rows
df_cut = df_cut[df_cut.start<>0]

After running:

>>> df_cut
     chr     start      stop geneID  tdiff
0  chr13  32889584  32889814  BRCA2    722
1  chr13  32890536  32890737  BRCA2   2457
2  chr13  32893194  32893466  BRCA2     -0

Answer 3

I modified the bed_prepare function to check for overlaps in previous and next genomic regions:

def bed_prepare(inp_bed):
    ''' Takes pandas dataframe bed file and identifies which regions overlap '''
    inp_bed['next_start'] = inp_bed['start'].shift(periods=-1)
    inp_bed['distance_to_next'] = inp_bed['next_start'] - inp_bed['stop']
    inp_bed['next_region_overlap'] = inp_bed['next_start'] <= inp_bed['stop']
    inp_bed['previous_stop'] = inp_bed['stop'].shift(periods=1)
    inp_bed['distance_from_previous'] = inp_bed['start'] - inp_bed['previous_stop']
    inp_bed['previous_region_overlap'] = inp_bed['previous_stop'] >= inp_bed['start']
    intermediate_bed = inp_bed
    return intermediate_bed

And then I used the Boolean outputs from these to do the variable storing for the writing step:

# Create empty dataframe to fill with parsed values                                                                   
new_bed = pd.DataFrame(data=np.zeros((0,len(columns))),columns=columns,dtype=int)

def bed_collapse(intermediate_bed, new_bed,columns=columns):
    ''' Takes a pandas dataframe bed file with overlap information and returns                                        
    genomic regions without overlaps '''
    output_row = []
    for row in bed.itertuples():
        output = {}
        if row[7] == False and row[10] == False:
            # If row doesn't overlap next row, insert into new dataframe unchanged.                                   
            output_row = list(row[1:5])
        elif row[7] == True and row[10] == False:
            # Only next region overlaps; take the chromosome and start coordinate                                     
            output_row = list(row[1:3])
        elif row[7] == True and row[10] == True:
            # Next and previous regions overlap. Skip row.                                                            
            pass
        elif row[7] == False and row[10]  == True:
            # Only previous region overlaps; append stop coordinate and geneID to output_row variable                 
            output_row.append(row[3])
            output_row.append(row[4])
        if row[7] == False:
            #Zip columns and output_row values together to form a dict for appending                                  
            for k, v in zip(columns,output_row): output[k] = v
            #print output                                                                                             
            new_bed = new_bed.append(output,ignore_index=True)
    output_bed = new_bed
    return output_bed

This has now solved my problem and gives the desired output specified in the question. :)

Answer 4

I will try to give you some pointers.

One pointer is that you want the get the rows based on a Series consisting of booleans that is shifted. Probably you can get a new shifted Series using:

Boolean_Series = intermediate_bed.loc[:,'next_region_overlap'].shift(periods=1, freq=None, axis=0, **kwds)

More background about this function: http://pandas.pydata.org/pandas-docs/dev/generated/pandas.DataFrame.shift.html

Second pointer is that by using this shifted Series you can get your Dataframe by:

int_bed = bed.loc[Boolean_Series, :] 

More about indexing can be found here: http://pandas.pydata.org/pandas-docs/dev/indexing.html

These are only pointers now, I do not know if this an actual working solution.