登陆服务器
终端——新建远程连接——用户：311lab——IP：10.2.12.189——连接进入输入密码：HELAB308——输入命令：cd /data1——输入：cd 311lab——输入：mkd wq , cd wq——mkd analysis mkd wq_rnaseq_file
文件位置：/data1/311lab/wq
去掉rawdata的双接头
在fastp软件所在的文件夹内
./fastp --in1 /data1/311lab/wq/wq_rnaseq_file/XXX0001_1.fq.gz --in2 /data1/311lab/wq/wq_rnaseq_file/XXX0001_2.fq.gz --out1 /data1/311lab/wq/analysis/XXX0001_1.fq.gz --out2 /data1/311lab/wq/analysis/XXX0001_2.fq.gz -l 20 -n 15 -q 5 -u 50 --json /data1/311lab/wq/analysis/0001.json --html /data1/311lab/wq/analysis/0001.html
hisat2比对index
index文件所在路径：/home/311lab/ref_V4_ensembl/index/hisat2/hisat
在analysis文件夹下新建文件夹：
mkd read_aln
双端测序（PE）
nohup hisat2 -x /home/311lab/ref_V4_ensembl/index/hisat2/hisat -1 /data1/311lab/wq/analysis/XXX0002_1.fq.gz -2 /data1/311lab/wq/analysis/XXX0002_2.fq.gz -S /data1/311lab/wq/analysis/read_aln/0002.sam（保存到路径文件夹下） &
samtools转化sam文件至bam文件
nohup Samtools view -bS XXXX.sam -o XXXX.bam &
文件排序：
nohup samtools sort XXX.bam -o XXX.srt.bam &
featurecounts计数
注释文件GTF所在路径：/home/311lab/ref_V4_ensembl/Zea_mays.B73_RefGen_v4.50.chr.gtf
nohup featureCounts -T 5 -p -t exon -g gene_id -a /home/311lab/ref_V4_ensembl/Zea_mays.B73_RefGen_v4.50.chr.gtf -o all.id.txt /data1/311lab/jim_file/practice/analysis/read_aln/srt.bam/*.bam &
FileZilla导入到本地电脑
主机IP：用户名：密码：端口：22
将all.id.txt文件传输到本地电脑上
igv可视化
topvlevel文件下载：ensembl玉米右侧栏——FASTA——DNA——dna.topleve——下载
导入toplevel文件查看每条染色体碱基，导入srt.bam文件看read比对的位置和密集程度，导入gtf文件看基因组遍布在染色体的位置。
R语言的下游分析画图
DESeq2的差异表达基因分析
数据的准备:
getwd()
setwd('/Users/wangqi/Documents/RNAseq/wqcounts')
count_tab = read.table('all.id.txt', header = T)
rownames(count_tab) = count_tab[,1]#行名改为第一列
count_tab = count_tab[, -c(1)]#去掉第一列
count_tab = (count_tab)[,-c(1:5)]#删掉前5列
colnames(count_tab)<- c("WT1","WT2","WT3","mutation1","mutation2","mutation3")#修改列名
colData <- read.table("samples.info.txt", header = T)#读取sample文件，取名为colData
str(counttab)
str(colData)
colData$sample=as.factor(colData$sample)#将colData数据框第一列数据修改为factor
colData$condition=as.factor(colData$condition)#将colData数据框第二列数据修改为factor
GEG_DESeq2分析:
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_tab, colData = colData, design= ~ condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_WT_vs_mutation")
#这个地方突变体较野生型相比的差异（后比前），不可以颠倒
res = res[order(res$padj),]#给res文件排序，用order（）命令，对行排序，对res的padj行排序
resDF = as.data.frame(res)#转换成可以看的数据框
write.csv(resDF, file = "wq_RNAseq_DEG.csv")
# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_WT_vs_mutation", type="apeglm")
PCA作图：
library(ggplot2)
vsd = vst(dds, blind = FALSE)
pcaData <- plotPCA(vsd, intgroup=c("condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
plotFig = ggplot(pcaData, aes(PC1, PC2, color=condition)) +#用代码保存PDF
#ggplot(pcaData, aes(PC1, PC2, color=condition)) ##不用代码直接输出
geom_point(size=3) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed()
ggsave(plotFig,filename = "wq_RNAseq_PCA.pdf")#代码保存PDF
MA-plot作图：
pdf("wq_RNAseq_MAplot.pdf")
plotMA(res, ylim=c(-2,2))
dev.off()
在Linux环境下分析差异表达
awk '$6 < 0.05 && $3 < -1' wq_RNAseq_DEG.txt 下调基因
awk '$6 < 0.05' wq_RNAseq_DEG.txt |wc -l查看有多少差异基因
awk '$6 < 0.05' wq_RNAseq_DEG.txt |less -S查看差异基因名称```
