| 1 | RNA=read.table('/virdir/Backup/run_1_gene_counts/combined_gene_count_run_1.txt',sep='\t',header=TRUE) |
|---|
| 2 | # save(RNA,file='/virdir/Backup/run_1_gene_counts/combined_gene_count_run_1_R_object') |
|---|
| 3 | # load('/virdir/Backup/run_1_gene_counts/combined_gene_count_run_1_R_object') |
|---|
| 4 | |
|---|
| 5 | ## get # zero counts per genes |
|---|
| 6 | |
|---|
| 7 | ind=which(RNA==0,arr.ind = TRUE) |
|---|
| 8 | h=hist(ind[,1],breaks=0:dim(RNA)[1],plot=FALSE) |
|---|
| 9 | nz=h$counts # for each gene, # samples with zero counts |
|---|
| 10 | |
|---|
| 11 | length(which(nz==0)) |
|---|
| 12 | [1] 3 |
|---|
| 13 | length(which(nz==1)) |
|---|
| 14 | [1] 326 |
|---|
| 15 | length(which(nz==2)) |
|---|
| 16 | [1] 15692 |
|---|
| 17 | |
|---|
| 18 | RNA_0=RNA[which(nz==2),] |
|---|
| 19 | |
|---|
| 20 | ## get IDs amsterdam sample |
|---|
| 21 | |
|---|
| 22 | id=as.character(read.table('/home/rjansen/NTR_ID_mapping.txt',sep='\t')[,2]) |
|---|
| 23 | m=match(id,colnames(RNA));length(which(is.na(m))) |
|---|
| 24 | |
|---|
| 25 | RNA_0_A = RNA_0[,m] # Amsterdam samples |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | ## spearman correlation |
|---|
| 29 | |
|---|
| 30 | Cs=cor(RNA_0,method='spearman') |
|---|
| 31 | mCs= apply(Cs, 2, median,na.rm = TRUE) |
|---|
| 32 | mCsa=mCs[m] # Amsterdam |
|---|
| 33 | mCsna=mCs[setdiff(1:length(mCs),m)] # non Amsterdam |
|---|
| 34 | |
|---|
| 35 | # Within a-dam: |
|---|
| 36 | Cs=cor(RNA_0_A,method='spearman') |
|---|
| 37 | mCsa_w = apply(Cs, 2, median) |
|---|
| 38 | |
|---|
| 39 | t.test(mCsna,mCsa) |
|---|
| 40 | t.test(mCsna,mCsa_w) |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | png('Histogram_Mean_correlations_A-dam_non_A-dam_Spear_0_within') |
|---|
| 44 | hist(mCsa_w,30, col="blue",freq=FALSE,xlim=c(0.8,1),xlab='Spearman correlations, for Amsterdam computed only within Amsterdam ',main='') |
|---|
| 45 | hist(mCsna,300, col='#FF000080' ,add=TRUE,transparant=1,freq=FALSE) |
|---|
| 46 | legend('top',c('A-dam','Non A-dam'),bty='n',fill=c("blue",'#FF000080')) |
|---|
| 47 | dev.off() |
|---|
| 48 | |
|---|
| 49 | |
|---|
| 50 | png('Histogram_Mean_correlations_A-dam_non_A-dam_Spear_0') |
|---|
| 51 | hist(mCsa,30, col="blue",freq=FALSE,xlim=c(0.8,1),xlab='Spearman correlations',main='') |
|---|
| 52 | hist(mCsna,300, col='#FF000080' ,add=TRUE,transparant=1,freq=FALSE) |
|---|
| 53 | legend('top',c('A-dam','Non A-dam'),bty='n',fill=c("blue",'#FF000080')) |
|---|
| 54 | dev.off() |
|---|
