Workflow(2) Differential expression

Posted on Edited on In

Workflow(2) Differential expression

Preparation

1.Download R, R studio

https://cran.r-project.org/mirrors.html
https://posit.co/download/rstudio-desktop/

Codes

1.Just copy the following codes and paste it into Rstudio

建议用txt等文本文件保存最原始的以及自己修改过的各种版本，方便改进
代码中所有需要更改的地方已经用”#”标注,”#”为R语言的注释符号
It will be better to save the original version and every version that you changed by your self to improve the reuse and accuracy of the code.
The place that can be changed has been annotated by “#”. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
      library("lattice")
      library(ggplot2)  
      library(corrplot)   
      setwd("F:/Taogroup/People/ZH/20210427_CVSP_WHOLE PROTEIN_PEAKS_18/result/")
      # 此处改为上一节最后存储下来的文件的文件路径 
      # This is the place to put the file path that got from former chapter
      DataHS <- read.table("proBeforeFCT.txt", header = T, sep = "\t")  
      #引号内填写自己的文件名称
      #Put file name in the " "
      rowN <- nrow(DataHS)  
      for(i in c(1 : rowN))  
      {  
        g1ave <- mean(as.numeric(DataHS[i, 1:3]))  
        g2ave <- mean(as.numeric(DataHS[i, 4:6]))  
        # 1：3，4：6分别表示两个组别的第一列数和最后一列数，根据自己的数据修改,g2ave2为后面一组
        #The number here represent the first and last column number in two different groups. Change it depends on yourself data
        FC <- g2ave - g1ave  
         if(FC > 1 || FC < -1)  
        {  
          DataHS[i, "fold.change"] <- FC  
          DataHS[i, "FC.significance"] <- "+"  
        }  
        else  
        {  
          DataHS[i, "fold.change"] <- FC  
          DataHS[i, "FC.significance"] <- ""  
        }  
      }  

      for(i in c(1 : rowN))  
      {  
       res <- t.test(as.numeric(DataHS[i, 1:3]), as.numeric(DataHS[i, 4:6]), alternative = "two.sided", conf.level = 0.9)
       #此处[]内数值需要更改为自己的列数
       #change the number in [] to yourselves data column number
       pVal <- res$p.value
       if(pVal < 0.05 || pVal > 0.95)
       {
         DataHS[i, "t.test.p.value"] <- -1*log10(pVal)
         DataHS[i, "t.test.significance"] <- "+"
       }
       else
       {
         DataHS[i, "t.test.p.value"] <- -1*log10(pVal)
         DataHS[i, "t.test.significance"] <- ""
       }
     }

     for(i in c(1 : rowN))
     {
       if(DataHS[i, "t.test.p.value"] > -1*log10(0.05)  && (DataHS[i, "fold.change"] >1 || DataHS[i, "fold.change"] < -1))
       {
         if(DataHS[i, "fold.change"] > 0)
         {
           DataHS[i, "sig"] = "Conup"
         }
         else
         {
      DataHS[i, "sig"] = "Tup"
    }
    #此处""内容根据组别更改，第一个为g2高表达，第二个为g1高表达  
    #"Conup" and "Tup" should be changed by your data, the first one represent the g2 up-regulated,the second one represent that g1 up-regulated
       }
       else
       {
         DataHS[i, "sig"] = "NoSignificance"
       }
     }
     write.table(DataHS, "T-FC-prodiff.txt", row.names = FALSE, quote = FALSE, sep = "\t")
     #引号内为差异表达的结果，可直接复制粘贴到excel中查看\  
     #"T-FC-prodiff.txt" is the file that contain all the information and data, you can check it in your workfile which was set at the begining.  
```     

### 2.Simple volcano plot and heatmap 
#### 2.1最简单的火山图可以由1.的数据直接绘制，所以这里介绍一下方法
首先运行以下两行  
Run the following two lines code after 1. 
```R

     A <- max(DataHS$fold.change)
     B <- min(DataHS$fold.change)

在Rstudio右侧查看A和B的数值，下面需要用到
Check the number in the right of Rstudio, it will be used later.
查看右侧AB数值 

1
2
3
4
5
6
7
8
9
10
11
12
ggplot(DataHS,aes(fold.change,t.test.p.value, color = sig))+
  geom_point()+xlim(-5,5)+labs(x="log2(FoldChange)", y="-log10(pValue)") + 
  scale_color_manual(values =c("red","gray","blue"))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)+
  theme_bw()+
  theme(legend.title = element_blank())
 dev.off()
 #第二行xlim后括号里的数值根据AB绝对值的最大值填写，这里是设置横坐标
 #第三行颜色部分建议修改，直接网上搜索R语言颜色代码，修改替代，另外出图时颜色顺序是按照分组的首字母排序的，有时候高低表达和不显著的颜色会交换，出现问题，建议在R里先看一下，有需要交换顺序的，看好要交换的颜色，将上述代码中的颜色位置直接换一下，再导出图片。导出图片点击右边plot放大，export，另外还可以修改代码直接导出pdf或者其他格式，具体方法会在下一节介绍。
 #The number in xlim(-5,5) is the length of x axis
 #The color in line3 could be changed, just search the color code for R, you will see many.

2.2 简易热图

将需要作图的数据复制粘贴出来为一个单独的txt文件，每一行为一个蛋白，每一列为一个样本，蛋白名放在最后一列 

1
2
3
4
5
6
data=read.table(file="proZ.txt",header=T,sep="\t")
rownames(data)=data[,9]
#[]内逗号后填写蛋白名的列数
pheatmap(data[,1:8],show_rownames =T,cluster_cols = F,fontsize_row=7)
#[]内逗号后填写蛋白表达量的列数，冒号表示从第几列到第几列，1：8表示从第一列到第八列
#同样，运行完这三行，会在Rstudio的右边出现图片，点击export导出

pheatmap功能很多，修改图片比例，注释字体，颜色，可以直接百度搜索”pheatmap参数” 学习
推荐：https://www.jianshu.com/p/31a92d6dc746