BOL: Owner

Calculate confidence interval for multiple columns in a matrix (table) using R function CI_normal and CI_tdist

EagleEye — Fri, 24 Jun 2016 18:01:04 -0500

#USAGE: CI_normal(matrix[,c(x:y)],97.5) - multiple column from a normal distribution
#USAGE: CI_tdist(matrix[,c(x:y)],97.5) - multiple column from a t distribution
#USAGE: CI_normal(matrix[,x],97.5) - single column from a normal distribution
#USAGE: CI_tdist(matrix[,x],97.5) - single column from a t distribution
#Created and updated by Santhilal Subhash on 2016/06/07 
CI_normal <- function(li,stat)
{
	cat(paste0("CI","\t","column","\t","Lower.limit","\t","Upper.limit","\n"))
	if(length(colnames(li))>1)
	{
		for(i in 1:length(colnames(li)))
		{
			sample <- colnames(li)
			error <- qnorm(stat/100)*sd(li[,i])/sqrt(length(li[,i]))
			left <- mean(li[,i])-error
			right <- mean(li[,i])+error
			cat(paste0(stat,"% norm dist","\t",sample[i],"\t",left,"\t",right,"\n"))

		}
	}
	else{
		error <- qnorm(stat/100)*sd(li)/sqrt(length(li))
		left <- mean(li)-error
		right <- mean(li)+error
		cat(paste0(stat,"% norm dist","\tNA","\t",left,"\t",right,"\n"))

	}

}
CI_tdist <- function(li,stat)
{
	cat(paste0("CI","\t","column","\t","Lower.limit","\t","Upper.limit","\n"))
	if(length(colnames(li))>1)
	{
		for(i in 1:length(colnames(li)))
		{
			sample <- colnames(li)
			error <- qt(stat/100,df=length(li[,i])-1)*sd(li[,i])/sqrt(length(li[,i]))
			left <- mean(li[,i])-error
			right <- mean(li[,i])+error
			cat(paste0(stat,"% t-dist","\t",sample[i],"\t",left,"\t",right,"\n"))

		}
	}
	else{

		error <- qt(stat/100,df=length(li)-1)*sd(li)/sqrt(length(li))
		left <- mean(li)-error
		right <- mean(li)+error
		cat(paste0(stat,"% t-dist","\tNA","\t",left,"\t",right,"\n"))

	}

}

RPKM normalization - Perl

EagleEye — Fri, 24 Jun 2016 17:53:09 -0500

#!/usr/bin/perl
#use strict;
#use warnings;
#Author: Santhilal Subhash
#Contact: santhilal.subhash@gu.se
#RPKM for RNAseq V1.3
#USAGE for sample input provided: perl rpkm_script_beta.pl sample_count_test.count 2:9 15 > sample_count_test.rpkm
#USAGE: perl rpkm_script_beta.pl input_count_file.txt ActualColumnStart:ActualColumnEnd ColumnGeneLength > results.rpkm
open $fh1, '<', $ARGV[0] or die $!;
open $fh2, '<', $ARGV[0] or die $!;
$total = 0;
$count = 0;
$cols=$ARGV[1];
$len_col=$ARGV[2];
while (<$fh1>) 
{
	@array=split("\t");
	($cstart,$cend)=split(":",$cols);
	for($i=$cstart-1;$i<=$cend-1;$i++)
	{
		$libarray[$i] += $array[$i];
	}
}

while (<$fh2>) 
{
	$next = <>;
	if ($next =~ /^#/) 
	{
		$header=$next;
		$header =~ s/#//;
		print "$header";
	}
	if ($next !~ /^#/) 
	{

 		@array=split("\t");
  		($cstart,$cend)=split(":",$cols);
  		for($i=$cstart-1;$i<=$cend-1;$i++)

		{
			if($array[$i]!=0)
			{
				$array_rpkm[$i]=((1000000000*$array[$i])/($libarray[$i]*$array[$len_col-1]));
			}
			else
			{
				$array_rpkm[$i]=0;
			}
		}
			local $" = "\t";
				
		print "$array[$cstart-2]@array_rpkm\t@array[$cend..$#array]";
	}
}

RPKM normalization- R script

EagleEye — Fri, 24 Jun 2016 17:45:55 -0500

P<-read.table("input_table.txt",sep="\t",header=T)

len <- ncol(p)

rownames(p) <- p[,1]

for(i in 2:ncol(p)-1)

{

d <-p[,i]

l <- p[,len] #accessing the length column

cS <- sum(as.numeric(p[,i])) #Total mapped reads per sample 

rpkm[[i]] <- (10^9)*(as.numeric(p[,i]))/(as.numeric(l)*cS)

rpkm[[1]] <- p[[1]]

}

write.table(rpkm,"output_table_rpkm.txt",sep="\t",quote=F,row.names=F)