linux基本命令
高通量相關(guān)基礎(chǔ)知識(shí)
1、 Anaconda
Anaconda是一個(gè)Python下和Canopy類(lèi)似的的科學(xué)計(jì)算環(huán)境，但用起來(lái)更加方便。自帶的包管理器conda也很強(qiáng)大，可以解決大多數(shù)軟件的安裝和配置
-
1-1、conda的安裝
Conda 官網(wǎng) https://conda.io/docs/index.html
mkdir bio_soft ## 創(chuàng)建一個(gè)文件夾，主演們用來(lái)放置軟件
cd bio_soft/ ## 進(jìn)入到該文件夾目錄下面
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 ## 將所下載的文件移動(dòng)到conda目錄下
cd conda
bash Anaconda2-5.1.0-Linux-x86_64.sh ## 運(yùn)行bash命令，一路回車(chē)或者yes，會(huì)自動(dòng)在用戶(hù)所屬環(huán)境創(chuàng)建一個(gè)anaconda2的文件夾，并且自動(dòng)添加環(huán)境變量
export PATH=/home/u883604/anaconda2/bin:$PATH ## 添加路徑到環(huán)境變量（已添加，可省略）
source ~/.bashrc
conda ##檢測(cè)是否安裝成功
添加幾個(gè)通道
conda config --add channels r
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
無(wú)論是conda默認(rèn)的軟件源還是bioconda軟件源都是國(guó)外的，速度非常慢，所以需要增加國(guó)內(nèi)軟件源，同時(shí)bioconda已經(jīng)有清華，中科大兩個(gè)國(guó)內(nèi)鏡像，也添加進(jìn)去
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --set show_channel_urls yes ## 設(shè)置搜索時(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的幾個(gè)常用命令
conda search PACKAGENAME ## 查找是否有所需的安裝文件
## 也可以在線(xiàn)查詢(xún) https://anaconda.org/
conda install PACKAGENAME ## 安裝所需軟件
## conda默認(rèn)安裝最新版本的軟件，若想安裝非最新版的軟件輸入：conda search bwa查看可選版本
conda install bwa=版本號(hào)
conda list ## 查看所有安裝的軟件
conda update 軟件名 可以對(duì)軟件進(jìn)行升級(jí):eg. conda update bwa
conda remove 卸載已經(jīng)安裝的軟件
更多詳細(xì)命令
-
1-3 、 Conda環(huán)境管理
## 由于我們下載的是python2.7版本的，但是如果我們需要運(yùn)行python3.0以上的呢？
# 創(chuàng)建一個(gè)名為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 ## 刪除一個(gè)已有的環(huán)境
參考鏈接
2、 RNA-seq軟件安裝
- Tophat2
conda 直接安裝
conda install -c bioconda sra-tools
conda install fastqc ## 不知道是網(wǎng)速還是怎么下載中斷好幾次，所以改為手動(dòng)安裝了
conda install trimmomatic
conda install tophat2
conda install bowtie2
conda install samtools
conda install cufflinks
有些軟件用conda 無(wú)法安裝，可以嘗試手動(dòng)安裝
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)境變量
## 臨時(shí)添加環(huán)境變量
export PATH=/home/u883604/bio_soft/tophat2/:$PATH
## 長(zhǎng)久使用需要修改 ~/.bashrc 文件
vim ~/.bashrc
PATH=/home/u883604/bio_soft/tophat2/:$PATH ## tophat2所在路徑
source ~/.bashrc
tophat2 -h
源代碼安裝
如果以上方法安裝不成功，可以使用源文件（文件名通常含src）的版本安裝
wget 下載鏈接
tar -zxvf 文件名
cd 文件
#看readme文件里的安裝提示，沒(méi)有該文件的話(huà)，自己網(wǎng)上查找
解壓命令：
tar -zxvf xxx.tar.gz
tar -jxvf xxx.tar.bz2
安裝分三步：
1. 配置： ./configure --prefix=想要指定的路徑
2. 編譯： make
3. 安裝： make install
若安裝出問(wèn)題(沒(méi)有權(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下載水稻的參考基因組文件和注釋文件
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ù)放置后臺(tái)運(yùn)行
下載擬南芥的參考基因組文件和注釋文件
## 創(chuàng)建文件夾用來(lái)放置參考基因組或注釋文件
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
下載測(cè)序數(shù)據(jù)
# 創(chuàng)建一個(gè)文件夾，用來(lái)存放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ù)上述返回?cái)?shù)字，獲取文件第33列，然后下載fastq文件
# tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}
# 也可以編寫(xiě)批量下載數(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
下載測(cè)序數(shù)據(jù)
# 創(chuàng)建一個(gè)文件夾，用來(lái)存放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)建一個(gè)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號(hào)/SRR號(hào).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)建一個(gè)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)建一個(gè)fastq-dump.sh文件，內(nèi)容如下
#PBS -N fastq-dump ## 指定任務(wù)名 （自定義）
#PBS -l nodes=1:ppn=4 ## 1表示一個(gè)線(xiàn)程 4表示4個(gè)核 同時(shí)也可以指定某個(gè)節(jié)點(diǎn)如 nodes=node32:ppn=4 （一般不需要改，4個(gè)基本夠）
#PBS <A8>Cl walltime=1:00:00 //請(qǐng)求任務(wù)執(zhí)行時(shí)間
#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
## 在這里我們是采用的是雙端測(cè)序的數(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
## 會(huì)得到8個(gè)fastq文件
得到了fastq文件我們就可以采用不同的RNA-seq protocol來(lái)進(jìn)行分析了
- hppRNA—a Snakemake-based handy parameter-free pipeline for RNA-Seq analysis of numerous samples
RNA-seq protocol compare.png
在進(jìn)行數(shù)據(jù)分析前必須進(jìn)行數(shù)據(jù)質(zhì)量的檢測(cè)，今天略過(guò)。。。。
Tophat –> Cufflink –> Cuffdiff
1、 對(duì)基因組序列建立索引
## 特別注意這里使用的index的gtf或者gff文件一定要與基因組序列文件要對(duì)應(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 ./ ##建立軟鏈接（類(lèi)似windows里面的快捷鍵）
2、將測(cè)序數(shù)據(jù)的reads比對(duì)到參考基因組上
mkdir align
cd align
## 這里使用的命令皆為首先創(chuàng)建sh文件，然后提交到集群運(yùn)行，具體格式參考之前
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ù)說(shuō)明
-p 8 表示使用8個(gè)線(xiàn)程
-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
## 注意如果是鏈特異性測(cè)序比對(duì)的時(shí)候需要加參數(shù) --library-type （后跟你的測(cè)序類(lèi)型）
fr-firststrand,fr-secondstrand
(詳情見(jiàn) http://m.itdecent.cn/p/a63595a41bed )
3、對(duì)每個(gè)樣本進(jìn)行轉(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)錄本進(jìn)行組裝
- 首先創(chuàng)建一個(gè)assemblies.txt（里面應(yīng)包含上一部得到的gtf文件的路徑）
## 創(chuàng)建文件可以使用vim assemblies.txt (將gft文件所在路徑粘貼進(jìn)去即可)
/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)錄本合并到一個(gè)轉(zhuǎn)錄本中，形成一個(gè)新的轉(zhuǎn)錄本注釋文件
## 不需要提交服務(wù)器 直接運(yùn)行即可
## 基因組注釋質(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ù)說(shuō)明
-g 后面跟參考基因組的注釋文件
-s 后面跟參考基因組序列文件
-p 線(xiàn)程數(shù)
5、獲得基因表達(dá)文件
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ù)說(shuō)明
grep -v 輸出不含inf的所有信息
Subread -> featureCounts -> DESeq2
1、 subread-buildindex 建立索引 （速度賊快 ~2.4min,可在計(jì)算節(jié)點(diǎn)直接運(yùn)行，不需提交）
subread-buildindex -o rice rice.fa
# rice.fa 為之前一樣的水稻參考基因組序列
# -o 輸出索引的前綴
2、 subjunc 進(jìn)行比對(duì)
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ù)說(shuō)明 (具體說(shuō)明直接運(yùn)行 subjunc )
-T 線(xiàn)程數(shù)
-i 之前建立的基因組索引的前綴
-r 接fastq文件
-R 雙端測(cè)序的fastq文件
-o 輸出名
subread-align.png
subread_align.png
3、 featureCounts 進(jìn)行定量
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 默認(rèn)exon
-g 默認(rèn)gene_id
用之前請(qǐng)查看自己的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ù)中只有第一列和從第七列(即第一個(gè)bam文件)到最后一列信息才是我們所需要的，即我們這里有四個(gè)bam文件所以提取第7到10列
# 參數(shù)說(shuō)明
\t 表示以制表符分割開(kāi)來(lái)
5、 將矩陣導(dǎo)入R中，采用DESeq2進(jìn)行差異分析
> countdata <- read.table("WT_osdrm2_matrix.out", header=TRUE, row.names=1) #導(dǎo)入數(shù)據(jù)
>
> head(countdata) # 查看數(shù)據(jù)前幾行(列名 太長(zhǎng)自己展示看)
>
> 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)建一個(gè)condition數(shù)據(jù)框
> dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition) ##構(gòu)建dds矩陣(后面很多分析都是基于這個(gè)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進(jìn)行繪制火山圖
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
原文獻(xiàn)圖
PP_原始火山圖.png
7、 GO富集分析
- 打開(kāi)網(wǎng)頁(yè)PCSD輸入我們得到的差異基因，這里拿上調(diào)的差異基因?yàn)槔?，看到有幾個(gè)選項(xiàng)，因?yàn)槲疫@里是沒(méi)有區(qū)分TE和Non TE的所以我直接選了第一個(gè)
GO分析—1.png
- 點(diǎn)擊submit
## 由于下載分GOSlim文件的ID為轉(zhuǎn)錄本ID，所以要先進(jìn)行處理
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，可能會(huì)報(bào)錯(cuò)，以下可能不適合。
ah <- AnnotationHub()
ah$species[which(ah$species == "Oryza sativa")] ## 水稻有幾個(gè)數(shù)據(jù)庫(kù)
query(ah, "Oryza sativa") ## 查找編號(hào)
> 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')
# 注意下面這個(gè)編號(hào)是會(huì)變的，所以你要輸入你自己上面查到的 `org.Oryza_sativa_Japonica_Group.eg.sqlite` 對(duì)應(yīng)的編號(hào)
rice_ref <- ah[['AH59059']]
## 查看rice_ref信息
str(rice_ref)
mode(rice_ref)
class(rice_ref)
columns(rice_ref) ## 這個(gè)可看
keytypes(rice_ref) ## 這個(gè)后面需要用到，這里我導(dǎo)入的是的GO號(hào)
head(keys(rice_ref,keytype = "GO"))
library("clusterProfiler")
PP_GO <- read.table(file.choose(),header=FALSE) ## 只有GO號(hào)的文件
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ù)圖
## 散點(diǎn)圖
## 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è)置熱圖每個(gè)格子的寬高
cluster_col = F,cluster_row= F, ## 按行還是按列聚類(lèi)，一般按行，即基因
main="The PP_data diff_gene heatmap ", (標(biāo)題)
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 ## 計(jì)算相關(guān)系數(shù)值
corrplot(cor_data,type = "upper") ## 繪制相關(guān)系數(shù)圖
## counts散點(diǎn)圖(也可以用來(lái)看相關(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 無(wú)重復(fù)樣本的差異分析
DEXseq
limma
GFOLD ## 無(wú)重復(fù)樣本的差異分析
最后我平時(shí)參考的論壇
github 一個(gè)開(kāi)源的網(wǎng)站，大佬們都把源碼放置在上面，比如已發(fā)表paper的源碼
感謝愛(ài)分享的各位生信人
逛的太多-- 就不一一舉例了
一直相信勤能補(bǔ)拙
天道酬勤.png
如有問(wèn)題請(qǐng)留言。。。
[https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.11](https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.11)
