|
10
|
1 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
|
|
|
2 args <- commandArgs(trailingOnly = TRUE)
|
|
|
3
|
|
|
4 # Read inputs
|
|
|
5 merged <- read.table( args[1] )
|
|
|
6 read_counts <- read.table( args[2] )
|
|
|
7
|
|
|
8 barcodes_nt <- merged[ do.call( paste, as.list( merged[ ,1:3 ] ) ) %in% do.call( paste, as.list( read_counts[ ( read_counts[ , 4 ] ) <= summary( read_counts[ , 4 ] )[ 5 ], 1:3 ] ) ) , 4 ]
|
|
|
9
|
|
|
10 ##make the matrix with the nucleotide freqs per position:
|
|
|
11 nt_counts <- matrix( nrow = 4, ncol = as.numeric(args[4]) )
|
|
|
12 for( h in 1:ncol( nt_counts ) ){
|
|
|
13 j <- 1
|
|
|
14 for( nt_local in c( "A","C","G","T" ) ) {
|
|
|
15 nt_counts[ j, h ] <- sum( substr( as.character( barcodes_nt ), h, h) == nt_local )
|
|
|
16 j <- j + 1
|
|
|
17 }
|
|
|
18 }
|
|
|
19
|
|
|
20 # Calculate frequencies
|
|
|
21 nt_freqs <- nt_counts / colSums( nt_counts )
|
|
|
22
|
|
|
23 nt_values <- list()
|
|
|
24 for( i in 1:ncol( nt_freqs ) ) {
|
|
|
25 nt_values[[ i ]] <- nt_freqs[ , i ]
|
|
|
26 }
|
|
|
27
|
|
|
28
|
|
|
29 all_posible_comb <- expand.grid( nt_values )
|
|
|
30
|
|
|
31 probs <- apply( all_posible_comb, 1, prod )
|
|
|
32
|
|
|
33 ###Create Mf_to_Sf:
|
|
|
34 results <- c()
|
|
|
35
|
|
|
36 i <- 1
|
|
|
37 results[ i ] <- sum( 1 - ( 1 - probs )**i )
|
|
|
38
|
|
|
39 while( results[ i ] <= as.numeric(args[3]) ) {
|
|
|
40 i <- i + 1
|
|
|
41 results[ i ] <- sum( 1 - ( 1 - probs )**i ) #Mf to Sf
|
|
|
42 }
|
|
|
43
|
|
|
44 #assign molecules number to unique barcode number:
|
|
|
45 Uf_to_Mf <- c()
|
|
|
46 for( i in 1:floor( max( results ) ) ) {
|
|
|
47 abs( results - i ) -> difference
|
|
|
48 Uf_to_Mf[ i ] <- which( difference == min( difference ) ) #if you want to know how many molecules underlie n unique barcodes ask Uf_to_Mf[n]
|
|
|
49 }
|
|
|
50
|
|
|
51 write( Uf_to_Mf, file = args[5] )
|
