Skip to main content

CAGE / RNA-seq normalizaton


Various normalization techniques in different tools

1. DESeq


Input a N by C matrix:
 N row ;  N genes
C column; Each column represents replica and different condition
Each cell : Represents the integer tag count value OR whatever value you want.

Algo:
1. From count matrix, for each row, calculate the geometric mean.
2. Diving each value by GM of corresponding row which generates size factor (SF) of each column.
3. Divide each count value by SF

MATLAB code

function doDESEQnorm

mymat=[   
    0    0    0    0    0    0    1;
    92    161    76    70    140    88    70;
    5    1    0    0    4    0    0;
    0    2    1    2    1    0    0 ];
   
% This shows four gene; in tow condition untreated and treated; First four replica for "untreated" and last %three column for "treated";


sizeRow = size(mymat,1);
sizeCol = size(mymat,2);

%% estimateSizeFactors %% 1 - Find Geometric Mean (GM) of Each row
gmRow = zeros(1,sizeRow);
for i=1:sizeRow
   
    curRow = mymat(i,:)';
    nzId = find(curRow);
    tmpVal = curRow(nzId);
    gm = geomean(tmpVal);
    gmRow(i) = gm;
   
end
mymat3= mymat;

%% estimateSizeFactors %% 2 - Divide each colum by corresponding row GM
afterDiv = mymat3;
for i=1:sizeRow
    afterDiv(i,:) = mymat3(i,:) / gmRow(i);
end



%% estimateSizeFactors %% 3 - SizeFactor: Take the median of Non-Zero Mormalized values
sizeFactor = zeros(1,sizeCol);
for col=1:sizeCol
   nzRow = find( afterDiv(:,col) );
   nzVal = afterDiv(nzRow,col);
   nzVal = sort(nzVal);
  
%    numberNZ = size(nzVal,1);
%    if rem(numberNZ, 2) ==0
%        idx = numberNZ/2 ;
%    else
%        idx = fix(numberNZ/ 2) +1;
%    end
%    sizeFactor(col)= nzVal(idx,1);

   sizeFactor(col) = median(nzVal);
  
end
sizeFactor

%% counts(sf, normalize=T) ; %% Do Normalzation
normValue = mymat;
for col=1:sizeCol
    normValue(:,col) = normValue(:,col)./sizeFactor(col);   
end

end

Comments

Popular posts from this blog

MATLAB cross validation

// use built-in function samplesize = size( matrix , 1); c = cvpartition(samplesize,  'kfold' , k); % return the indexes on each fold ///// output in matlab console K-fold cross validation partition              N: 10    NumTestSets: 4      TrainSize: 8  7  7  8       TestSize: 2  3  3  2 ////////////////////// for i=1 : k    trainIdxs = find(training(c,i) ); %training(c,i);  // 1 means in train , 0 means in test    testInxs  = find(test(c,i)       ); % test(c,i);       // 1 means in test , 0 means in train    trainMatrix = matrix (  matrix(trainIdxs ), : );    testMatrix  = matrix (  matrix(testIdxs  ), : ); end //// now calculate performance %%  calculate performance of a partition     selectedKfoldSen=[];selectedKfoldSpe=[];selectedKfoldAcc=[];     indexSen=1;indexSpe=1;indexAcc=1;     if ( kfold == (P+N) )% leave one out         sensitivity = sum(cvtp) /( sum(cvtp) + sum(cvfn) )         specificity = sum(cvtn) /( sum(cvfp) + sum(cvtn) )         acc

R tutorial

Install R in linux ============ In CRAN home page, the latest version is not available. So, in fedora, Open the terminal yum list R  --> To check the latest available version of r yum install R --> install R version yum update R --> update current version to latest one 0 find help ============ ?exact topic name (  i.e.   ?mean ) 0.0 INSTALL 3rd party package  ==================== install.packages('mvtnorm' , dependencies = TRUE , lib='/home/alamt/myRlibrary/')   #  install new package BED file parsing (Always use read.delim it is the best) library(MASS) #library(ggplot2) dirRoot="D:/research/F5shortRNA/TestRIKEN/Rscripts/" dirData="D:/research/F5shortRNA/TestRIKEN/" setwd(dirRoot) getwd() myBed="test.bed" fnmBed=paste(dirData, myBed, sep="") # ccdsHh19.bed   tmp.bed ## Read bed use read.delim - it is the  best mybed=read.delim(fnmBed, header = FALSE, sep = "\t", quote = &q

MATLAB confusion matrix

%  test_class  & predicted_class must be same dimension % 'order' - describes the order of label. Here labels are 'g' as positive and 'h' as negative [C,order] = confusionmat( test_class(1: noSampleTest), predicted_class, 'order', ['g' ;'h'] ) tp = C(1,1); fn = C(1,2); fp = C(2,1); tn = C(2,2); sensitivity = tp /( tp + fn ) specificity = tn /( fp + tn ) accuracy = (tp+tn) / (tp+fn+fp+tn) tpr = sensitivity fpr = 1-specificity precision = tp /( tp + fp ) fVal = (2*tpr*precision)/(tpr+precision)