From ad2e4bd07f8b94f9a3c058e98efd6ec05ad3e7cb Mon Sep 17 00:00:00 2001 From: Monika Cechova Date: Fri, 8 Mar 2024 13:34:24 -0800 Subject: [PATCH] Fixes bug at chromosome boundaries, refactors the code, het_ratio is now the average across all second and first alleles --- HetDetection.R | 167 +++++++++++++++++++++---------------------------- 1 file changed, 72 insertions(+), 95 deletions(-) diff --git a/HetDetection.R b/HetDetection.R index 178cf29..bac2dd7 100644 --- a/HetDetection.R +++ b/HetDetection.R @@ -1,105 +1,82 @@ -library(ggplot2) -library(directlabels) -library(plyr) -require(gridExtra) -require(scales) -library(RColorBrewer) -library(reshape2) -library(data.table) -library(dplyr) -library(tidyr) -library(zoo) -​ -###change the working directory (2nd line) and saved file name (last line) before running### -​ -#clear variables -rm(list=ls(all=TRUE)) -​ -#set working directory -setwd("~/Documents/Eichler Lab/Weekly plans/T2Tv1/Hets") -​ +#This script filters the output by nucfreq and reports "regions where the second most common base was present in at least 10% of reads in at least 5 positions within a 500 bp region" + +require(tidyr) +require(data.table) +require(dplyr) + +# Check if any command line arguments are provided +if(length(commandArgs(trailingOnly = TRUE)) == 0) { + stop("No filename provided. Please provide a filename as a command-line argument.") +} + +# Get the filename from the command line +filename <- commandArgs(trailingOnly = TRUE)[1] + + #load NucFreq bed file -df = fread(file.choose(), stringsAsFactors = FALSE, fill=TRUE, quote="", header=FALSE, skip=2) +df = read.table(filename, stringsAsFactors = FALSE, quote="", header=TRUE) cols = c("chr", "start", "end", "first", "second") colnames(df) <- cols -​ + #determine the ratio of the first and second most common bases df$het_ratio = round(df$second/(df$first+df$second)*100, 1) -​ + #filter if the het ratio is >= 10% -df1 = df %>% +df = df %>% group_by(chr) %>% filter(het_ratio >= 10) -​ -#calculate the distance (in bp) between consecutive positions -df2 = df1 %>% - group_by(chr) %>% - mutate(distance = start - lag(start, default = start[1])) -​ -#shift the distance column up one row -shift <- function(x, n){ - c(x[-(seq(n))], rep(NA, n)) -} -​ -df2$distance <- shift(df2$distance, 1) -​ -#filter rows with a distance <=500 bp between positions (i.e. the het must have another base change within 500 bp) -df3 = df2 %>% - group_by(chr) %>% - filter(distance <= 500) -df3 = df3 %>% - group_by(chr) %>% - mutate(distance2 = start - lag(start, default = start[1])) -df3$distance2 <- shift(df3$distance2, 1) -​ -#duplicate top row and change its distance to 501 bp (to get rows in register) -df4 = df3 %>% - group_by(chr) %>% - filter(row_number() <= 1) %>% - bind_rows(df3) -df5 = df4 %>% - arrange(start, .by_group = TRUE) %>% - mutate(distance2 = replace(distance2, row_number() == 1, 501)) -​ -#shift up the end column to get the range of the hets on one row -df5$end <- shift(df5$end, 1) -​ -#filter only if there are 5 consecutive rows of distance <=500 bp (i.e. the het must have 5 base changes within 500 bp) -r <- with(with(df5, rle(distance2<=500)),rep(lengths,lengths)) -df5$het <- with(df5,distance2<=500) & (r>=4) -​ -#filter the row if it contains a het -df6 <- filter(df5, het == "TRUE") -​ -###determine the max and min coordinates of a region with a distance <= 500 bp -​ -#get the min and max coordinates -het_df = df6 %>% - group_by(chr) %>% - mutate(distance3 = start - lag(start, default = start[1])) -het_df$distance3 <- shift(het_df$distance3, 1) -het_df2 = het_df %>% - group_by(chr) %>% - filter((distance3 >= 500) | lag(distance3 >= 500) | (row_number() >= (n())) | (row_number() == 1) ) -​ -#shift the end column to have the min and max coordinates on one row -het_df2$end <- shift(het_df2$end, 1) -​ -#take every other row (1, 3, etc.) -het_df3 = het_df2 %>% - group_by(chr) %>% - filter(row_number() %% 2 == 1) -​ -#take the differences of the min and max coordinates -het_df3 = het_df3 %>% + +#calculate the distances betwen variants +df <- df %>% + arrange(chr, start) %>% group_by(chr) %>% - mutate(het_length = end - start) -​ -#remove those with negative lengths -het_df_filtered = het_df3 %>% - select(chr, start, end, het_ratio, het_length) %>% - filter(het_length > 0) -​ -#print the table -write.table(het_df_filtered, "all.hets.winnowmap.greaterThan10_2.tbl", row.names = F, quote = F, sep="\t") + mutate( + lag_distance = lead(start) - end + ) + +#is distance smaller than 500 bp? if yes, set to true +df$closby<-df$lag_distance<=500 + +#the last rows are going to have NA, since there are more variants than distances +#let's fill it with the value of the previous row +df <- df %>% + fill(closby, .direction = "down") + +#count how many times you get consecutive TRUE or FALSE (looking for the consecutive small distances) +#this is achieved using run-length encoding, and then expanding it so that the number of rows matches +df$rle<-rep(rle(df$closby)$lengths,rle(df$closby)$lengths) + +#assign a unique group id for each cluster, since the groups should not get mixed up +df$group<-rleid(paste0(df$closby,df$rle)) + +#only keep rows that belong to clusters +df<-df[!df$closby==FALSE,] + +#only keep rows that cluster at least 5 variants +df<-df[df$rle>=5,] + +het_ratio_per_group <- df %>% + group_by(group) %>% + summarise(sum_second = sum(second), + sum_first = sum(first)) %>% + mutate(het_ratio = round(sum_second/(sum_first+sum_second)*100, 1)) + + +df_final <- df %>% + group_by(chr,group) %>% + filter(row_number() == 1 | row_number() == n()) %>% + summarise(start = first(start), + end = last(end)) + +#calculate the lengths of het regions +df_final$het_length<-df_final$end-df_final$start + +#merge with the het frequency information +df_table<-merge(df_final,het_ratio_per_group,by="group") + +#only keep start and end columns +df_table<-df_table[,c("chr","start","end","het_ratio","het_length")] +df_table<-na.omit(df_table) +write.table(df_table, paste0(filename,"all.hets.filtering.tbl"), row.names = F, quote = F, sep="\t") +