Introduction

Let’s QC this data.

library("gplots")
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
library("reshape2")

Read

tmp <- read.table("3col.tsv")
x <- t(acast(tmp, V1~V2, value.var="V3"))

# chr length
chr <- read.table("pc_chr.tsv")

Contigs analysis

head(x[order(-rowMeans(x)),],20)
##                         SR2       SR5       SR6       SR9       UR3      UR4
## JAFJYM010000023.1 60034.100 65067.200 71149.500 30582.100 67633.600 79037.90
## JAFJYM010000037.1 32391.100 33789.900 36596.600 35418.100 35145.800 51067.90
## JAFJYM010000077.1 11014.800 10608.300 11512.000 11391.500 10788.500 15041.60
## JAFJYM010000001.1  3636.710  4424.240  4886.950  6411.380  4885.240  7458.48
## JAFJYM010000123.1  1155.410  1874.400  1816.260  2664.090  2637.950  2568.20
## JAFJYM010000002.1  1537.120  1711.730  1940.950  2375.790  1915.320  2438.88
## JAFJYM010000056.1  1429.570  1330.410  1525.020  3070.230  1467.130  2284.33
## JAFJYM010000046.1  1359.860  1042.370  1304.300  2666.030  1244.940  1873.89
## JAFJYM010000133.1  1380.570  1367.860  1618.590  1708.720  1317.300  2085.19
## JAFJYM010000015.1  1001.370  1315.320  1472.330  1044.590  1166.920  1514.83
## JAFJYM010000132.1   931.604  1037.090  1182.180   784.173  1354.330  1997.88
## JAFJYM010000083.1   942.801   942.805   975.445  1054.970  1269.890  1452.49
## JAFJYM010000064.1   749.488  2065.430  1002.360  1012.630   842.926  1639.51
## JAFJYM010000016.1   989.119  1259.020  1088.760  1040.450  1026.830  1326.80
## JAFJYM010000050.1   982.120   686.993   672.868   793.630   759.910  1181.09
## JAFJYM010000054.1   551.335   843.594   669.686   664.807   642.681  1199.92
## JAFJYM010000027.1   514.032   725.846   757.607   729.130   484.618  1028.65
## JAFJYM010000084.1   548.178   753.210   598.329   625.742   535.672  1097.17
## JAFJYM010000100.1   272.615   417.597  1054.100   562.873   429.181  1164.50
## JAFJYM010000022.1   338.669   506.854   693.319   771.554   412.870   800.29
##                         UR7       UR8
## JAFJYM010000023.1 64005.100 41281.800
## JAFJYM010000037.1 41711.200 39977.200
## JAFJYM010000077.1 11026.400 12357.600
## JAFJYM010000001.1  3967.430  5419.440
## JAFJYM010000123.1  2484.550  2297.500
## JAFJYM010000002.1  2049.730  2235.890
## JAFJYM010000056.1  1463.240  2514.140
## JAFJYM010000046.1  1405.640  1978.980
## JAFJYM010000133.1  1358.350  1800.820
## JAFJYM010000015.1   789.036  1577.650
## JAFJYM010000132.1   852.391   904.777
## JAFJYM010000083.1  1295.980  1095.080
## JAFJYM010000064.1   505.647  1105.370
## JAFJYM010000016.1  1022.760  1099.770
## JAFJYM010000050.1   866.358  1281.710
## JAFJYM010000054.1   485.482   658.868
## JAFJYM010000027.1   456.110   736.670
## JAFJYM010000084.1   434.933   566.731
## JAFJYM010000100.1   282.846   894.592
## JAFJYM010000022.1   374.637   751.081
chr <- chr[order(chr$V1),]

xn <- x/chr$V2 * 1000000

Differential analysis

res <- apply(X=xn,MARGIN=1,function(i){
  g1 <- i[1:4]
  g2 <- i[5:8]
  xt <- t.test(x=g1,y=g2,alternative="two.sided")
  xp <- xt$p.value
  xlfc <- unname(log2(xt$estimate[2]/xt$estimate[1]))
  return(c(xp,xlfc))
})

res <- t(res)

colnames(res) <- c("pval","lfc")

res <- as.data.frame(res)

res$fdr <- p.adjust(res$pval)

top <- head(res[order(res$pval),],20)

top
##                          pval        lfc       fdr
## JAFJYM010000052.1 0.003959024  1.0001092 0.4909190
## JAFJYM010000055.1 0.005819693  0.7735410 0.7158223
## JAFJYM010000083.1 0.020489087  0.3849057 1.0000000
## JAFJYM010000094.1 0.036305039  0.8365210 1.0000000
## JAFJYM010000082.1 0.044161051  0.5205446 1.0000000
## JAFJYM010000048.1 0.073855101  1.4151266 1.0000000
## JAFJYM010000072.1 0.075622379  1.4343284 1.0000000
## JAFJYM010000087.1 0.076379745  1.8073549 1.0000000
## JAFJYM010000025.1 0.080535590  0.8763054 1.0000000
## JAFJYM010000026.1 0.090549255  0.4905615 1.0000000
## JAFJYM010000004.1 0.096990115  0.4425570 1.0000000
## JAFJYM010000037.1 0.109091063  0.2809075 1.0000000
## JAFJYM010000031.1 0.135970645 -1.1724087 1.0000000
## JAFJYM010000123.1 0.136608726  0.4113811 1.0000000
## JAFJYM010000112.1 0.143313675  1.9444387 1.0000000
## JAFJYM010000058.1 0.157825693  0.5826446 1.0000000
## JAFJYM010000050.1 0.160031601  0.3830255 1.0000000
## JAFJYM010000073.1 0.175001121  0.6513053 1.0000000
## JAFJYM010000099.1 0.208508123  1.2017802 1.0000000
## JAFJYM010000104.1 0.213486048  2.3258879 1.0000000
for (i in 1:6) { 
  r=which(rownames(xn) == rownames(top)[i])
  barplot(xn[r,],main=rownames(xn)[r])
}

Session info

sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] reshape2_1.4.4 gplots_3.1.1  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7         gtools_3.9.2       digest_0.6.27      bitops_1.0-7      
##  [5] plyr_1.8.6         R6_2.5.0           jsonlite_1.7.2     magrittr_2.0.1    
##  [9] evaluate_0.14      highr_0.9          KernSmooth_2.23-20 rlang_0.4.11      
## [13] stringi_1.7.3      jquerylib_0.1.4    bslib_0.2.5.1      rmarkdown_2.10    
## [17] tools_4.1.1        stringr_1.4.0      xfun_0.25          yaml_2.2.1        
## [21] compiler_4.1.1     caTools_1.18.2     htmltools_0.5.1.1  knitr_1.33        
## [25] sass_0.4.0