RNA-seq完整分析過程

linux基本命令

高通量相關(guān)基礎(chǔ)知識

1骗爆、 Anaconda

Anaconda是一個Python下和Canopy類似的的科學(xué)計算環(huán)境次氨，但用起來更加方便。自帶的包管理器conda也很強大摘投，可以解決大多數(shù)軟件的安裝和配置
1-1煮寡、conda的安裝


Conda 官網(wǎng) https://conda.io/docs/index.html

mkdir bio_soft ## 創(chuàng)建一個文件夾，主演們用來放置軟件

cd bio_soft/   ## 進入到該文件夾目錄下面

wget https://repo.continuum.io/archive/Anaconda2-5.1.0-Linux-x86_64.sh ## 下載python2.7的Anaconda 

mkdir conda

mv Anaconda2-5.1.0-Linux-x86_64.sh conda ## 將所下載的文件移動到conda目錄下

cd conda

bash Anaconda2-5.1.0-Linux-x86_64.sh ## 運行bash命令犀呼，一路回車或者yes幸撕，會自動在用戶所屬環(huán)境創(chuàng)建一個anaconda2的文件夾锥咸，并且自動添加環(huán)境變量

export PATH=/home/u883604/anaconda2/bin:$PATH ## 添加路徑到環(huán)境變量（已添加租悄，可省略）

source ~/.bashrc 

conda ##檢測是否安裝成功




添加幾個通道

conda config --add channels r

conda config --add channels defaults

conda config --add channels conda-forge

conda config --add channels bioconda




無論是conda默認的軟件源還是bioconda軟件源都是國外的掸读，速度非常慢淘太，所以需要增加國內(nèi)軟件源，同時bioconda已經(jīng)有清華馏锡，中科大兩個國內(nèi)鏡像示括，也添加進去

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/

conda config --set show_channel_urls yes  ## 設(shè)置搜索時顯示通道地址

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/msys2/

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/

 查看目前conda軟件源情況

conda info

1-2览效、 Conda的幾個常用命令


conda search PACKAGENAME ## 查找是否有所需的安裝文件

## 也可以在線查詢 https://anaconda.org/

conda install PACKAGENAME  ## 安裝所需軟件

## conda默認安裝最新版本的軟件罪佳，若想安裝非最新版的軟件輸入：conda search bwa查看可選版本 

conda install bwa=版本號

conda list ## 查看所有安裝的軟件

conda update 軟件名 可以對軟件進行升級:eg. conda update bwa

conda remove 卸載已經(jīng)安裝的軟件

更多詳細命令

conda commands

1-3 逛漫、 Conda環(huán)境管理


## 由于我們下載的是python2.7版本的，但是如果我們需要運行python3.0以上的呢赘艳？

# 創(chuàng)建一個名為python3的環(huán)境酌毡，指定Python版本是3.6

conda create --name python3 python=3.6

source activate python3 ## 加載python3的環(huán)境

source deactivate ## 使用結(jié)束克握，可以退出環(huán)境

conda remove --name python34 --all ## 刪除一個已有的環(huán)境

參考鏈接

生信札記

PeterYuan

2、 RNA-seq軟件安裝

Tophat2

conda 直接安裝


conda install -c bioconda sra-tools

conda install fastqc  ## 不知道是網(wǎng)速還是怎么下載中斷好幾次阔馋，所以改為手動安裝了

conda install trimmomatic

conda install tophat2

conda install bowtie2

conda install samtools

conda install cufflinks

有些軟件用conda 無法安裝玛荞，可以嘗試手動安裝

tophat.png

預(yù)編譯版本安裝


cd bio_soft

[tophat2官網(wǎng)](http://ccb.jhu.edu/software/tophat/index.shtml)

wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.1.1.Linux_x86_64.tar.gz ## 下載已經(jīng)編譯好了的

mv tophat-2.1.1.Linux_x86_64/ tophat2

cd tophat2/

## 添加環(huán)境變量

## 臨時添加環(huán)境變量

export PATH=/home/u883604/bio_soft/tophat2/:$PATH 

## 長久使用需要修改 ~/.bashrc 文件

vim ~/.bashrc

    PATH=/home/u883604/bio_soft/tophat2/:$PATH ## tophat2所在路徑

source ~/.bashrc   

tophat2 -h

源代碼安裝

如果以上方法安裝不成功，可以使用源文件（文件名通常含src）的版本安裝


wget 下載鏈接

tar -zxvf 文件名

cd 文件

#看readme文件里的安裝提示呕寝，沒有該文件的話，自己網(wǎng)上查找




解壓命令：

tar -zxvf xxx.tar.gz

tar -jxvf xxx.tar.bz2

安裝分三步：

1. 配置： ./configure --prefix=想要指定的路徑

2. 編譯： make

3. 安裝： make install

若安裝出問題(沒有權(quán)限在系統(tǒng)目錄下安裝)婴梧，需要清空上次編譯的文檔

make clean

./configure --prefix=安裝路徑

make 

make insatll #可執(zhí)行文件通常在bin目錄下下梢，環(huán)境變量設(shè)置和源碼安裝一樣

FastQC安裝


wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip

unzip fastqc_v0.11.7.zip

cd FastQC/

chmod 777 fastqc

./fastqc -h

vim ~/.bashrc 

    export PATH=/home/u883604/bio_soft/FastQC/:$PATH

source ~/.bashrc

cufflinks安裝


wget http://cole-trapnell-lab.github.io/cufflinks/assets/downloads/cufflinks-2.2.1.Linux_x86_64.tar.gz

tar -xzvf cufflinks-2.2.1.Linux_x86_64.tar.gz

cd cufflinks-2.2.1.Linux_x86_64/

vim ~/.bashrc

    export PATH=/home/u883604/bio_soft/cufflinks-2.2.1.Linux_x86_64/:$PATH

source ~/.bashrc

3、 RNA-seq分析流程

3.1下載水稻的參考基因組文件和注釋文件

水稻數(shù)據(jù)庫


mkdir -p rna_practice/ref 創(chuàng)建

nohup wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.con & ## 參考基因組

nohup wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.gff3 & ## 注釋文件

# 使用nohup &可以將任務(wù)放置后臺運行

下載擬南芥的參考基因組文件和注釋文件


## 創(chuàng)建文件夾用來放置參考基因組或注釋文件

mkdir -p rna_practice/data

cd rna_practice/ref

nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz &

nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa.gz &

nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gff3.gz &

nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gtf.gz &

# 解壓縮

gzip -d *.gz

下載測序數(shù)據(jù)


# 創(chuàng)建一個文件夾塞蹭，用來存放FASTQ文件(或者公司返回的數(shù)據(jù))

mkdir -p rna_practice/data/

cd rna_practice/data/




# 獲取樣本信息

wget http://www.ebi.ac.uk/arrayexpress/files/E-MTAB-5130/E-MTAB-5130.sdrf.txt

tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 32,36 |sort -u




# 下載數(shù)據(jù)

# fastq文件下載鏈接在第幾列

# head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI

# 根據(jù)上述返回數(shù)字孽江，獲取文件第33列，然后下載fastq文件

# tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}

# 也可以編寫批量下載數(shù)據(jù)的shell腳本番电，如下：

perl -alne 'if(/.*(ftp:.*gz).*/){print "nohup wget $1 &"}' E-MTAB-5130.sdrf.txt >fq_data_download.sh

bash fq_data_download.sh

paper.png

下載測序數(shù)據(jù)


# 創(chuàng)建一個文件夾岗屏，用來存放FASTQ文件(或者公司返回的數(shù)據(jù))

mkdir -p rna_practice/data/

cd rna_practice/data/




#第一種方法

nohup prefetch -v SRR1999345 &

nohup prefetch -v SRR1999346 &

nohup prefetch -v SRR1999347 &

nohup prefetch -v SRR1999348 &




#第二種方法

cat SRR_Acc_List.txt | xargs prefetch -v




#第三種方法

### 創(chuàng)建一個down.sh文件，里面內(nèi)容如下：

    #漱办！/bin/bash

    for i in `seq 5 8`

    do

    prefetch -v SRR199934${i} &

bash down.sh

本次數(shù)據(jù)


## ftp格式為SRR/SRR+SRR數(shù)字的前三位/SRR號/SRR號.sra

nohup wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503149/SRR3503149.sra &

nohup wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503154/SRR3503154.sra &

nohup wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503151/SRR3503151.sra &

nohup wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503152/SRR3503152.sra &

將sra文件轉(zhuǎn)換為fq文件


第一種：創(chuàng)建一個fastq-dump.sh文件这刷，內(nèi)容如下

    #!/bin/bash




    for i in `seq 5 9`

    do 

    nohup fastq-dump --split-3 --gzip SRR199934${i}.sra

    done&

nohup bash fastq-dump.sh &




第二種：創(chuàng)建一個fastq-dump.sh文件，內(nèi)容如下

#PBS -N fastq-dump   ## 指定任務(wù)名 （自定義）

#PBS -l nodes=1:ppn=4  ## 1表示一個線程  4表示4個核  同時也可以指定某個節(jié)點如 nodes=node32:ppn=4 （一般不需要改娩井，4個基本夠）

#PBS <A8>Cl walltime=1:00:00  //請求任務(wù)執(zhí)行時間

#PBS -q batch 

#PBS -V  

cd $PBS_O_WORKDIR

ls *.sra | xargs fastq-dump --split-3 --gzip  ## 所執(zhí)行的命令（自定義）




qsub fastq-dump.sh  ##提交任務(wù)到服務(wù)器




集群使用命令參考：

http://hpc.ncpgr.cn:8093/mediawiki/index.php/PBS%E4%BD%BF%E7%94%A8

qsub.png

qstat.png

本流程將sra文件轉(zhuǎn)換為fastq文件使用命令


fastq-dump 一般使用參數(shù)

## 單端數(shù)據(jù)使用命令

fastq-dump  SRR***.sra -O out_path

fastq-dump  --fasta  SRR306998.sra  -O out_path

## 雙端數(shù)據(jù)使用命令

fastq-dump --split-3  SRR***.sra -O out_path




## 在這里我們是采用的是雙端測序的數(shù)據(jù)暇屋，所以使用雙端的命令

fastq-dump --split-3 --gzip Cr-DJ-2-1.sra  # --gzip 是指生成的文件為壓縮的文件，可以節(jié)省占用內(nèi)存

fastq-dump --split-3 --gzip Cr-DJ-2-2.sra

fastq-dump --split-3 --gzip osdrm2-1.sra

fastq-dump --split-3 --gzip osdrm2-2.sra

## 會得到8個fastq文件

得到了fastq文件我們就可以采用不同的RNA-seq protocol來進行分析了

hppRNA—a Snakemake-based handy parameter-free pipeline for RNA-Seq analysis of numerous samples

RNA-seq protocol compare.png

在進行數(shù)據(jù)分析前必須進行數(shù)據(jù)質(zhì)量的檢測洞辣，今天略過咐刨。。扬霜。定鸟。

Tophat –> Cufflink –> Cuffdiff

1、對基因組序列建立索引


## 特別注意這里使用的index的gtf或者gff文件一定要與基因組序列文件要對應(yīng)

mkdir  rice_bowtie2_index  ## 存放水稻參考基因組文件

cd rice_bowtie2_index

bowtie2-build all.fa rice  ## rice為輸出的索引的前綴

mkdir rice_gtf  ## 存放水稻的注釋文件即gff3文件或者gtf文件

mkdir data 

cd data

ln -s /public/home/qliu/data/B512_data/RNA_data/osdrm2/*.fastq ./ ##建立軟鏈接（類似windows里面的快捷鍵）

2著瓶、將測序數(shù)據(jù)的reads比對到參考基因組上


mkdir align

cd align

## 這里使用的命令皆為首先創(chuàng)建sh文件联予，然后提交到集群運行，具體格式參考之前

tophat -p 8 -o CrDJ-1 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 /public/home/qliu/RNA_practice/rice_bowtie2_index/rice /public/home/qliu/RNA_practice/data/Cr-DJ-2-1_1.fastq /public/home/qliu/RNA_practice/data/Cr-DJ-2-1_2.fastq




tophat -p 8 -o CrDJ-2 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 /public/home/qliu/RNA_practice/rice_bowtie2_index/rice /public/home/qliu/RNA_practice/data/Cr-DJ-2-2_1.fastq /public/home/qliu/RNA_practice/data/Cr-DJ-2-2_2.fastq




tophat -p 8 -o osdrm-1 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 /public/home/qliu/RNA_practice/rice_bowtie2_index/rice /public/home/qliu/RNA_practice/data/osdrm2-1_1.fastq /public/home/qliu/RNA_practice/data/osdrm2-1_2.fastq




tophat -p 8 -o osdrm2-2 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 /public/home/qliu/RNA_practice/rice_bowtie2_index/rice /public/home/qliu/RNA_practice/data/osdrm2-2_1.fastq /public/home/qliu/RNA_practice/data/osdrm2-2_2.fastq




## 參數(shù)說明

-p 8 表示使用8個線程

-o CrDJ-1 表示結(jié)果文件輸出到CrDJ-1文件夾中（不需要自己提前創(chuàng)建）

-G 后面跟基因組注釋文件

/public/home/qliu/RNA_practice/rice_bowtie2_index/rice  表示之前所建立的基因組index所指定的前綴

所以參考的模板為

tophat -p 8 -G genes.gtf -o C1_R1_thout index C1_R1_1.fq C1_R1_2.fq




## 注意如果是鏈特異性測序比對的時候需要加參數(shù) --library-type （后跟你的測序類型）

fr-firststrand,fr-secondstrand 

(詳情見 http://www.reibang.com/p/a63595a41bed )

3蟹但、對每個樣本進行轉(zhuǎn)錄本組裝


mkdir cufflinks

cufflinks -p 8 -o CrDJ-1 /public/home/qliu/RNA_practice/align/CrDJ-1/accepted_hits.bam




cufflinks -p 8 -o CrDJ-2 /public/home/qliu/RNA_practice/align/CrDJ-2/accepted_hits.bam




cufflinks -p 8 -o osdrm2-1 /public/home/qliu/RNA_practice/align/osdrm-1/accepted_hits.bam




cufflinks -p 8 -o osdrm2-2 /public/home/qliu/RNA_practice/align/osdrm2-2/accepted_hits.bam

4躯泰、將所有的轉(zhuǎn)錄本進行組裝

首先創(chuàng)建一個assemblies.txt（里面應(yīng)包含上一部得到的gtf文件的路徑）


## 創(chuàng)建文件可以使用vim assemblies.txt (將gft文件所在路徑粘貼進去即可)

/public/home/qliu/RNA_practice/cufflinks/CrDJ-1/transcripts.gtf

/public/home/qliu/RNA_practice/cufflinks/CrDJ-2/transcripts.gtf

/public/home/qliu/RNA_practice/cufflinks/osdrm2-1/transcripts.gtf

/public/home/qliu/RNA_practice/cufflinks/osdrm2-1/transcripts.gtf

將所有的轉(zhuǎn)錄本合并到一個轉(zhuǎn)錄本中，形成一個新的轉(zhuǎn)錄本注釋文件


## 不需要提交服務(wù)器 直接運行即可
## 基因組注釋質(zhì)量高华糖，如果只拿差異基因麦向，不需要跑這一步
mkdir cuffmerge

cd cuffmerge

cuffmerge -g /public/home/qliu/RNA_practice/rice_gtf/all.gff3 -s /public/home/qliu/RNA_practice/rice_bowtie2_index/all.fa -p 8 assemblies.txt




參數(shù)說明

-g 后面跟參考基因組的注釋文件

-s 后面跟參考基因組序列文件

-p 線程數(shù)

5、獲得基因表達文件


mkdir cuffdiff

cuffdiff -o ./diff_out/ -b /public/home/qliu/RNA_practice/rice_bowtie2_index/all.fa -p 8 -L CrDJ,osdrm2 -u /public/home/qliu/RNA_practice/cuffmerge/merged_asm/merged.gtf /public/home/qliu/RNA_practice/align/CrDJ-1/accepted_hits.bam,/public/home/qliu/RNA_practice/align/CrDJ-2/accepted_hits.bam /public/home/qliu/RNA_practice/align/osdrm-1/accepted_hits.bam,/public/home/qliu/RNA_practice/align/osdrm2-2/accepted_hits.bam

6客叉、根據(jù)文章設(shè)立的閾值诵竭，挑選差異基因


awk '{if(($10<-2)&&($11<0.001))print $3"\t"$8"\t"$9"\t"$10}' gene_exp.diff | grep -v 'inf' > down.txt  ## 篩選出下調(diào)的基因（log2_fold_change < -2 & pvalue < 0.001）




awk '{if(($10>2)&&($11<0.001))print $3"\t"$8"\t"$9"\t"$10}' gene_exp.diff | grep -v 'inf' > up.txt ## 篩選出上調(diào)的基因（log2_fold_change > 2 & pvalue < 0.001




## 參數(shù)說明 

grep -v 輸出不含inf的所有信息

Subread -> featureCounts -> DESeq2

1话告、 subread-buildindex 建立索引（速度賊快 ~2.4min,可在計算節(jié)點直接運行，不需提交）


subread-buildindex -o rice rice.fa  

# rice.fa 為之前一樣的水稻參考基因組序列

# -o 輸出索引的前綴

2卵慰、 subjunc 進行比對


mkdkir align

cd align

subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice -r /public/home/qliu/RNA_practice/data/Cr-DJ-2-1_1.fastq -R /public/home/qliu/RNA_practice/data/Cr-DJ-2-1_2.fastq -o Cr-DJ_1 




subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice -r /public/home/qliu/RNA_practice/data/Cr-DJ-2-2_1.fastq -R /public/home/qliu/RNA_practice/data/Cr-DJ-2-2_2.fastq -o Cr-DJ_2




subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice -r /public/home/qliu/RNA_practice/data/osdrm2-1_1.fastq -R /public/home/qliu/RNA_practice/data/osdrm2-1_2.fastq -o osdrm2_1




subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice -r /public/home/qliu/RNA_practice/data/osdrm2-2_1.fastq -R /public/home/qliu/RNA_practice/data/osdrm2-2_2.fastq -o osdrm2_2




## 查看bam文件命令,比如查看前幾行

samtools view file.bam |head




## 參數(shù)說明 (具體說明直接運行 subjunc )

-T 線程數(shù)

-i 之前建立的基因組索引的前綴

-r 接fastq文件

-R 雙端測序的fastq文件

-o 輸出名

subread-align.png

subread_align.png

3沙郭、 featureCounts 進行定量


mkdir featureCounts

cd featureCounts

featureCounts -p -T 6 -a /public/home/qliu/RNA_practice/rice_gtf/rice.gtf -o Cr_DJ-osdrm2_fCount.out /public/home/qliu/RNA_practice/subread/align/Cr-DJ_1 /public/home/qliu/RNA_practice/subread/align/Cr-DJ_2 /public/home/qliu/RNA_practice/subread/align/osdrm2_1 /public/home/qliu/RNA_practice/subread/align/osdrm2_2 

## 注意參數(shù)

-t 默認exon

-g 默認gene_id

用之前請查看自己的gff文件或者gtf文件

4、提取定量的信息


awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' Cr_DJ-osdrm2_fCount.out > WT_osdrm2_matrix.out




## Cr_DJ-osdrm2_fCount.out數(shù)據(jù)中只有第一列和從第七列(即第一個bam文件)到最后一列信息才是我們所需要的裳朋，即我們這里有四個bam文件所以提取第7到10列




# 參數(shù)說明

\t 表示以制表符分割開來

5病线、將矩陣導(dǎo)入R中，采用DESeq2進行差異分析


> countdata <- read.table("WT_osdrm2_matrix.out", header=TRUE, row.names=1) #導(dǎo)入數(shù)據(jù)

> 

> head(countdata) # 查看數(shù)據(jù)前幾行(列名 太長自己展示看)

> 

> colnames(countdata) <- c("WT_1","WT_2","DRM2_1","DRN2_2") # 修改列名




> head(countdata)

               WT_1 WT_2 DRM2_1 DRN2_2

LOC_Os01g01010 1250 1250    948    574

LOC_Os01g01019   66   66    136    114

LOC_Os01g01030  535  535    665    464

LOC_Os01g01040 1531 1531   1225    847

LOC_Os01g01050  556  556    742    524

LOC_Os01g01060 2925 2925   3312   2374




> dim(countdata) # 查看數(shù)據(jù)維度鲤嫡，即幾行幾列

[1] 55986     4




>  condition <- factor(c("WT","WT","DRM2","DRM2")) # 賦值因子送挑，即變量




> library(DESeq2) # 這一步以后基本可以復(fù)制粘貼




> coldata <- data.frame(row.names=colnames(countdata), condition) # 創(chuàng)建一個condition數(shù)據(jù)框




> dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition) ##構(gòu)建dds矩陣(后面很多分析都是基于這個dds矩陣)




> dds <- DESeq(dds)




> resdata <- results(dds)




> table(resdata$padj<0.05) # p<0.05的基因數(shù)







> res_padj <- resdata[order(resdata$padj), ]  ##按照padj列的值排序




> names(resdata)[1] <- "Gene" # 將第一列的列名改為Gene




> write.table(resdata, file="diffexpr_padj_results.csv",sep = "\t",row.names = F) ## 將結(jié)果文件保存到本地




## 篩選差異基因

> subset(resdata,pvalue < 0.001) -> diff ## 先篩選pvalue < 0.01的行

 

> subset(diff,log2FoldChange < -2) -> down ## 篩選出下調(diào)的




> subset(diff,log2FoldChange > 2) -> up ## 篩選出上調(diào)的




> dim(down)

[1] 1124   11




> dim(up)

[1] 3023   11

6、 DESeq2分析中得到的resdata進行繪制火山圖


rm(list=ls())

## 加載DESeq2中生成的resdata文件

resdata <- read.csv(file.choose(),header = T , sep = "\t")

threshold <- as.factor(ifelse(resdata$padj < 0.001 & 

                                abs(resdata$log2FoldChange) >= 2 ,

                              ifelse(resdata$log2FoldChange >= 2 ,'Up','Down'),'Not'))

ggplot(resdata,aes(x=log2FoldChange,y=-log10(padj),colour=threshold)) +

  xlab("log2(Fold Change)")+ylab("-log10(qvalue)") +

  geom_point(size = 0.5,alpha=1) +

  ylim(0,200) + xlim(-12,12) +

  scale_color_manual(values=c("green","grey", "red"))

火山圖.png

原文獻圖

PP_原始火山圖.png

7暖眼、 GO富集分析

打開網(wǎng)頁PCSD輸入我們得到的差異基因惕耕，這里拿上調(diào)的差異基因為例，看到有幾個選項诫肠，因為我這里是沒有區(qū)分TE和Non TE的所以我直接選了第一個

GO分析—1.png

點擊submit


## 由于下載分GOSlim文件的ID為轉(zhuǎn)錄本ID司澎，所以要先進行處理

sed -e 's/\.[0-9]\{1\}//g' all.GOSlim_assignment > new_GOSlim.txt


## 將gene ID轉(zhuǎn)換為GO ID

all.GOSlim  <- read.table("/public/home/qliu/data/rice_data/rice_all.dir/new_GOSlim.txt",sep = "\t",header =T)




colnames(all.GOSlim) <- c("ID","GO_number","Background","Pathyway","IEA","TAIR")




up_gene <- read.table("up.gene",header=FALSE)




colnames(up_gene) <- "ID"




library(dplyr)




## 篩選出up gene 所在的GO_ID

left_join(up_gene,all.GOSlim,by = "ID") -> temp.data




GO_ID <- as.matrix(temp.data[,2])




GO_ID <- as.character(GO_ID[,1])




require(AnnotationHub)
# 在此，此包版本為 2.16.1栋豫，如果您 R 版本為 4.0挤安，可能會報錯，以下可能不適合笼才。


ah <- AnnotationHub()

ah$species[which(ah$species == "Oryza sativa")] ## 水稻有幾個數(shù)據(jù)庫

query(ah, "Oryza sativa") ## 查找編號

> query(ah, "Oryza sativa")

    AnnotationHub with 2 records

    # snapshotDate(): 2017-10-27 

    # $dataprovider: Inparanoid8, ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/

    # $species: Oryza sativa, Oryza sativa_Japonica_Group

    # $rdataclass: Inparanoid8Db, OrgDb

    # additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer,

    #   rdatadateadded, preparerclass, tags, rdatapath, sourceurl, sourcetype 

    # retrieve records with, e.g., 'object[["AH10561"]]' 

    title                                    

    AH10561 | hom.Oryza_sativa.inp8.sqlite             

    AH59059 | org.Oryza_sativa_Japonica_Group.eg.sqlite




## subset(ah, species == 'Oryza sativa' & rdataclass == 'OrgDb')
# 注意下面這個編號是會變的漱受，所以你要輸入你自己上面查到的  `org.Oryza_sativa_Japonica_Group.eg.sqlite` 對應(yīng)的編號
rice_ref <- ah[['AH59059']]







## 查看rice_ref信息

str(rice_ref) 

mode(rice_ref)

class(rice_ref)

columns(rice_ref)  ## 這個可看

keytypes(rice_ref) ## 這個后面需要用到，這里我導(dǎo)入的是的GO號

head(keys(rice_ref,keytype = "GO"))




library("clusterProfiler")

PP_GO <- read.table(file.choose(),header=FALSE)  ## 只有GO號的文件

PP_GO <- as.character(PP_GO[,1])

PP_test <- enrichGO(gene         = PP_GO,

                    OrgDb         = rice_ref,

                    keyType = "GO",

                    pAdjustMethod = "BH", ## “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”

                    ont      = "CC" , ## 可選"BP" "CC" "MF"

                    pvalueCutoff  = 0.01,

                    qvalueCutoff  = 0.05)

dotplot(PP_test, showCategory=30)

enrichMap(PP_test, vertex.label.cex=1.2,layout=igraph::layout.kamada.kawai)

barplot(PP_test,showCategory=12,font.size=7,title="groupGO Cellular Component")

library(topGO)

plotGOgraph(PP_test)

GO—1.png

GO-2.png

GO-3.png

GO-4.png

其他一些圖


## 熱圖  

## 相關(guān)系數(shù)圖

## 散點圖




## heatmap

rm(list=ls()) ## 清除當(dāng)前環(huán)境變量

resdata <- read.csv(file.choose(),header = T , sep = "\t") ## 導(dǎo)入數(shù)據(jù)PP_data_diffexpr_padj_results

subset(resdata,padj < 0.001 & abs(resdata$log2FoldChange) >= 2) -> diff_gene ## 篩選出差異基因

diff_gene[,8:ncol(diff_gene)] -> heatmap_data ## 提取counts value所在列

## ncol(diff_gene) ## 表示該數(shù)據(jù)有多少列

rownames(heatmap_data) <- as.array(diff_gene[,1]) ## 添加行名

as.matrix(heatmap_data) -> x  ## 轉(zhuǎn)換數(shù)據(jù)格式




library(pheatmap)  # 加載包 如果未安裝先install.packages("pheatmap")

pheatmap(x,scale="row",cellwidth=40,cellheight=0.1, ## 設(shè)置熱圖每個格子的寬高

         cluster_col = F,cluster_row= F, ## 按行還是按列聚類骡送，一般按行昂羡，即基因

         main="The PP_data diff_gene heatmap ", (標題)

         fontsize=5,treeheight_row = 2,show_rownames= F, ## 是否顯示行名(基因名)

         cutree_row=1,display_numbers = FALSE, ## 是否顯示數(shù)值

         color = colorRampPalette(c("green","white","red"))(100), ## 設(shè)置顏色

         clustering_distance_rows = "correlation", #filename="out.pdf"







## 相關(guān)系數(shù)圖

## 還是用上面的resdata數(shù)據(jù)

install.packages("corrplot")

library(corrplot)

resdata[,8:ncol(diff_gene)] -> heatmap_data  ## 提取WT以及osdrm2 counts value所在的列

cor(heatmap_data) -> cor_data ## 計算相關(guān)系數(shù)值

corrplot(cor_data,type = "upper") ## 繪制相關(guān)系數(shù)圖




## counts散點圖(也可以用來看相關(guān)性)

install.packages("ggplot2")

library(ggplot2)

ggplot(resdata,aes(x= resdata$WT_1,y=resdata$WT_2)) + 

  geom_point() + xlim(0,80000) + ylim(0,80000) +

  geom_abline() + theme_bw() +

  xlab("WT_1") +ylab("WT_2") +

  theme(panel.border = element_blank(),panel.grid.major = element_blank(),

        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))




ggplot(resdata,aes(x= resdata$DRM2_1,y=resdata$DRN2_2)) + 

  geom_point() + xlim(0,80000) + ylim(0,80000) +

  geom_abline() + theme_bw() +

  xlab("DRM2_1") + ylab("DRM2_2") +

  theme(panel.border = element_blank(),panel.grid.major = element_blank(),

        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))

heatmap.png

corrplot.png

WT.png

drm2.png

其他流程

HISAT2 -> HTSeq -> DESEq2
HISAT2 ->StringTie -> Ballgown
kallisto -> sleuth
....

差異分析R包

DESeq / DESeq2
edger
DEGSeq 無重復(fù)樣本的差異分析
DEXseq
limma
GFOLD ## 無重復(fù)樣本的差異分析

最后我平時參考的論壇

github 一個開源的網(wǎng)站，大佬們都把源碼放置在上面摔踱，比如已發(fā)表paper的源碼
生信技能樹
omicshare
rabbit gao個人博客
簡書
基因課
小木蟲
r-bloggers
感謝愛分享的各位生信人
逛的太多-- 就不一一舉例了

一直相信勤能補拙

天道酬勤.png

如有問題請留言虐先。。派敷。
[https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.11](https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.11)

最后編輯于：2020.07.17 15:34:05

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末蛹批，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子篮愉，更是在濱河造成了極大的恐慌腐芍，老刑警劉巖，帶你破解...
沈念sama閱讀 206,839評論 6贊 482
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件试躏，死亡現(xiàn)場離奇詭異猪勇，居然都是意外死亡，警方通過查閱死者的電腦和手機颠蕴，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,543評論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門泣刹，熙熙樓的掌柜王于貴愁眉苦臉地迎上來助析，“玉大人，你說我怎么就攤上這事椅您⊥饧剑” “怎么了？”我有些...
開封第一講書人閱讀 153,116評論 0贊 344
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵掀泳，是天一觀的道長雪隧。經(jīng)常有香客問我，道長开伏，這世上最難降的妖魔是什么膀跌？我笑而不...
開封第一講書人閱讀 55,371評論 1贊 279
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮固灵，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘劫流。我一直安慰自己巫玻，他們只是感情好，可當(dāng)我...
茶點故事閱讀 64,384評論 5贊 374
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布祠汇。她就那樣靜靜地躺著仍秤，像睡著了一般。火紅的嫁衣襯著肌膚如雪可很。梳的紋絲不亂的頭發(fā)上诗力，一...
開封第一講書人閱讀 49,111評論 1贊 285
城市分裂傳說
那天，我揣著相機與錄音我抠，去河邊找鬼苇本。笑死，一個胖子當(dāng)著我的面吹牛菜拓，可吹牛的內(nèi)容都是我干的瓣窄。我是一名探鬼主播，決...
沈念sama閱讀 38,416評論 3贊 400
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼纳鼎，長吁一口氣：“原來是場噩夢啊……” “哼俺夕！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起贱鄙，我...
開封第一講書人閱讀 37,053評論 0贊 259
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤劝贸，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后逗宁，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體映九，經(jīng)...
沈念sama閱讀 43,558評論 1贊 300
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 36,007評論 2贊 325
?白月光啟示錄
正文我和宋清朗相戀三年疙剑，在試婚紗的時候發(fā)現(xiàn)自己被綠了氯迂。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片践叠。...
茶點故事閱讀 38,117評論 1贊 334
活死人
序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖嚼蚀，靈堂內(nèi)的尸體忽然破棺而出禁灼，到底是詐尸還是另有隱情，我是刑警寧澤轿曙，帶...
沈念sama閱讀 33,756評論 4贊 324
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布弄捕，位于F島的核電站，受9級特大地震影響导帝，放射性物質(zhì)發(fā)生泄漏守谓。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點故事閱讀 39,324評論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一您单、第九天我趴在偏房一處隱蔽的房頂上張望斋荞。院中可真熱鬧，春花似錦虐秦、人聲如沸平酿。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,315評論 0贊 19
一樁弒父案悦陋，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽蜈彼。三九已至，卻和暖如春俺驶，著一層夾襖步出監(jiān)牢的瞬間幸逆，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 31,539評論 1贊 262
情欲美人皮
我被黑心中介騙來泰國打工暮现，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留还绘，地道東北人。一個月前我還...
沈念sama閱讀 45,578評論 2贊 355
代替公主和親
正文我出身青樓送矩，卻偏偏與公主長得像蚕甥，于是被迫代替她去往敵國和親。傳聞我的和親對象是個殘疾皇子栋荸，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 42,877評論 2贊 345