Normalisation de la couverture de profondeur parmi les échantillons

Question

Il s'agit d'une question ouverte qui vise à définir des états pour chaque position dans le génome (correspondant à un "CpG" sites) qui varient entre les échantillons. La raison de cette question est que les outils disponibles, définissent le statut d'une fenêtre CpG et non pour des CpG individuels. Les colonnes sont les suivantes: Nombre chromosomique (chr), position initiale (début) et finale (fin) d'une base d'intérêt, la couverture attendue (profondeur), la couverture observée à différents 6 animaux (profondeur1-profondeur6).

data <- "chr start end depth depth1 depth2 depth3 depth4 depth5 depth6 
chr1 3273 3273 7 200 35 1 200 850 0 
chr1 3274 3274 3 50 25 5 300 1500 2 
chr1 3275 3275 8 600 15 8 100 300 5 
chr1 3276 3276 4 30 2 10 59 20 0 
chr1 3277 3277 25 20 7 4 600 45 0" 
data <- read.table(text=data, header=T) 

J'ai besoin de définir une colonne avec les états de chaque ligne, les états sont: non couverts région, tour à tour couvert et souvent couvert. Pour ce faire, d'abord, j'ai besoin de faire une normalisation de la profondeur entre les échantillons pour obtenir des valeurs qui peuvent être comparées entre les individus. et, en second lieu, je dois définir la gamme entre les états (maintenant, toute gamme est acceptable);

J'ai trouvé ce script qui fait quelque chose de similaire à ce dont j'ai besoin pour normaliser les profondeurs, mais je ne pouvais pas encore l'appliquer à mon cas (ce script a été conçu pour "CpG windows" fonction ". dans chaque fenêtre

setMethod("normalizeCoverage", "methylRawList", 

         function(obj,method){ 

          if(method=="median"){ 
          x=sapply(obj,function(x) median(x$coverage)) 
          }else if(method=="mean"){ 
          x=sapply(obj,function(x) mean(x$coverage)) 
          }else{ 
          stop("method option should be either 'mean' or 'median'\n") 
          } 
          sc.fac=max(x)/x #get scaling factor 
          for(i in 1:length(obj)){ 
          all.cov=obj[[i]]$coverage 
          fCs =obj[[i]]$numCs/all.cov 
          fTs =obj[[i]]$numT/all.cov 
          obj[[i]]$coverage=round(sc.fac[i]*obj[[i]]$coverage) 
          obj[[i]]$numCs =round(obj[[i]]$coverage*fCs) 
          obj[[i]]$numTs =round(obj[[i]]$coverage*fTs) 
          } 
          obj 

    }) 

J'ai vérifié aussi ce « paquet bord » R, qui est utilisé pour les données de normalisation RNA-seq, qui ressemble à ceci:

calcNormFactors(object, method=c("TMM","RLE","upperquartile","none"), refColumn = NULL, 
     logratioTrim = .3, sumTrim = 0.05, doWeighting=TRUE, Acutoff=-1e10, p=0.75) 

mais aussi je ne pouvais pas appliquer à mes données encore

Ce que je souhaite pour mon résultat final est quelque chose comme ceci:

chr start State 
chr1 3273 Often 
chr1 3274 alternatively 
chr1 3275 no 
chr1 3276 often 
chr1 3277 no 

mais je serais vraiment satisfait que de la profondeur normalisée à chaque couverture de l'échantillon.

Answer 1

Pour la première partie du problème (normalisation)

« Calculer les valeurs de couverture ajustée profondeur, une application linéaire jusqu'à la fonction cpm à la matrice de comptage pourrait être suffisant. Cela va convertir vos comptes en nombre -par millions de valeurs que vous pouvez ensuite comparer de manière informelle entre les échantillons ". (Aaron Lun, Cambridge, Royaume-Uni)

Normaliser par "cpm" par paquet "déligneuse" de R

datamatrix <- data [, c(4:10)] 
library (edgeR) 

#grouping factor 
group <- c(1, 2, 2, 2, 2, 2, 2) #groups of each sample) 

#create a DGEList 
y <- DGEList(counts=datamatrix,group=group) 

######### 
$counts 
    depth depth1 depth2 depth3 depth4 depth5 depth6 
1  7 200  35  1 200 850  0 
2  3  50  25  5 300 1500  2 
3  8 600  15  8 100 300  5 
4  4  30  2  10  59  20  0 
5 25  20  7  4 600  45  0 

$samples 
     group lib.size norm.factors 
depth  1  47   1 
depth1  2  900   1 
depth2  2  84   1 
depth3  2  28   1 
depth4  2  1259   1 
depth5  2  2715   1 
depth6  2  7   1 
################## 
#normalize 
y <- calcNormFactors(y) 

######## 
> y 
An object of class "DGEList" 
$counts 
    depth depth1 depth2 depth3 depth4 depth5 depth6 
1  7 200  35  1 200 850  0 
2  3  50  25  5 300 1500  2 
3  8 600  15  8 100 300  5 
4  4  30  2  10  59  20  0 
5 25  20  7  4 600  45  0 

$samples 
     group lib.size norm.factors 
depth  1  47 0.7423735 
depth1  2  900 0.5526927 
depth2  2  84 0.9534847 
depth3  2  28 0.8652676 
depth4  2  1259 0.6588115 
depth5  2  2715 1.0358307 
depth6  2  7 4.3289213 
########################################## 

> cpm(y) 
     depth  depth1 depth2 depth3 depth4  depth5 depth6 
1 200621.61 402071.90 436993.56 41275.42 241125.49 302245.841  0.00 
2 85980.69 100517.97 312138.26 206377.10 361688.24 533375.014 66001.27 
3 229281.84 1206215.69 187282.96 330203.36 120562.75 106675.003 165003.16 
4 114640.92 60310.78 24971.06 412754.21 71132.02 7111.667  0.00 
5 716505.76 40207.19 87398.71 165101.68 723376.48 16001.250  0.00 

Normalisée!

Même avec la normalisation, j'ai 3 échantillons qui ont beaucoup de valeurs égales à zéro, car ils ont une faible couverture. Je pense que je vais devoir les supprimer de l'analyse. J'ai pensé faire un test PCA pour voir comment ces échantillons sont groupés.

Je voudrais des commentaires sur la méthode utilisée pour la normalisation et pour la deuxième partie de mon problème