R/FGI/rmd3/doplot.r

   1 #!/usr/bin/env littler
   2
   3 # Usage: Rscript gene_id_plot.R sample_report.csv
   4
   5 # 读入单个样本的报告表
   6 #cmargs = commandArgs(trailingOnly = T)
   7 sample_report_df = read.csv(argv[1])
   8 #sample_report_df = read.csv("C:\\Users\\wangjingdong\\Desktop\\gene_id_rscript\\input_csv_report.csv")
   9
  10 # 准备数据框为目标画图形式
  11 sample_report_df$End = sample_report_df$Pos
  12 for (rsid in rownames(sample_report_df)){
  13   genotype = sample_report_df[rsid,"GT"]
  14   if (genotype == "Null"){sample_report_df[rsid, "Val"] = 0}
  15   else{
  16     allele_1 = substr(genotype, 1, 1)
  17     allele_2 = substr(genotype, 2, 2)
  18     if (allele_1 == allele_2){sample_report_df[rsid, "Val"] = 1}
  19     else{sample_report_df[rsid, "Val"] = 0.5}
  20   }
  21 }
  22
  23 # 分成4个碱基文件
  24 base_group = function(df, base){
  25   bool = c()
  26   for (i in 1:length(rownames(df))){
  27     geno = df[i,2]
  28     if (geno == "Null"){bool = append(bool,FALSE)}
  29     else{
  30       alleles = c(substr(geno, 1, 1), substr(geno, 2, 2))
  31       if (base %in% alleles){bool = append(bool,TRUE)}
  32       else{bool = append(bool,FALSE)}
  33     }
  34   }
  35   return(df[rownames(df)[bool],c("Chr","Pos","End","Val")])
  36 }
  37
  38 matrix_A = base_group(sample_report_df, "A")
  39 matrix_A$PlotColor = "#5050FF"
  40
  41 matrix_T = base_group(sample_report_df, "T")
  42 matrix_T$PlotColor = "#CC9900"
  43
  44 matrix_G = base_group(sample_report_df, "G")
  45 matrix_G$PlotColor = "#00C000"
  46
  47 matrix_C = base_group(sample_report_df, "C")
  48 matrix_C$PlotColor = "#E00000"
  49
  50 #out.file <- "nCircos.pdf"
  51 out.file <- argv[2]
  52 #pdf(file=out.file, height=8, width=8)
  53 png(file=out.file, height=2160, width=2160,pointsize=48)
  54
  55 library(RCircos)
  56 data(UCSC.HG19.Human.CytoBandIdeogram)
  57 RCircos.Set.Core.Components(cyto.info=UCSC.HG19.Human.CytoBandIdeogram, tracks.inside=4, tracks.outside=0, chr.exclude=c("chrX", "chrY"))
  58
  59 rcircos.params <- RCircos.Get.Plot.Parameters()
  60 rcircos.params$track.background <- 'white'
  61 rcircos.params$grid.line.color <- 'white'
  62 rcircos.params$hist.width <- 160
  63 rcircos.params$track.padding <- 0
  64 RCircos.Reset.Plot.Parameters(rcircos.params)
  65
  66 RCircos.Set.Plot.Area(margins=0)
  67 legend("center", legend=c("A", "C", "G", "T"), col=c("#5050FF", "#E00000","#00C000","#CC9900"), lwd=16)
  68
  69 RCircos.Chromosome.Ideogram.Plot()
  70
  71 RCircos.Histogram.Plot(hist.data=matrix_A, data.col=4,track.num=1, side="in",is.sorted=F)
  72 RCircos.Histogram.Plot(hist.data=matrix_C, data.col=4,track.num=2, side="in",is.sorted=F)
  73 RCircos.Histogram.Plot(hist.data=matrix_G, data.col=4,track.num=3, side="in",is.sorted=F)
  74 RCircos.Histogram.Plot(hist.data=matrix_T, data.col=4,track.num=4, side="in",is.sorted=F)
  75
  76 dev.off()