###############################################RNA-seq pipline###############################################

##############################################准备工作#########################################################
#已安装软件路径：/lustre/home/acct-lctest/stu411/RNA-seq/softwares
drwxrwxr-x  9 stu411 stu411 16384 6月  13 08:26 bbmap
drwxrwxr-x  8 stu411 stu411  4096 6月  13 23:21 FastQC
drwxrwxr-x  5 stu411 stu411 12288 6月  13 23:45 hisat2-2.0.0-beta
drwxrwxr-x 11 stu411 stu411 12288 6月  14 00:06 samtools-1.10
drwxr-xr-x  6 stu411 stu411  4096 5月  13 11:46 subread-2.0.1-Linux-x86_64
这些软件的运行命令已经加入环境变量，直接运行即可

#原始数据
原始数据路径：/lustre/home/acct-lctest/stu411/RNA-seq/raw_data
数据说明：小鼠不同发育时间点细胞的RNA-seq双端测序数据。
time0，time1分别表示不同发育时间点，每个时间点分别有两个样品，表示为replicate1和replicate2，共4个样品。
1.9M -rw-r--r-- 1 stu411 stu411 1.9M 6月  14 00:29 time0_rep1_R1.fastq.gz
2.2M -rw-r--r-- 1 stu411 stu411 2.2M 6月  14 00:29 time0_rep1_R2.fastq.gz
2.0M -rw-r--r-- 1 stu411 stu411 2.0M 6月  14 00:29 time0_rep2_R1.fastq.gz
2.2M -rw-r--r-- 1 stu411 stu411 2.2M 6月  14 00:29 time0_rep2_R2.fastq.gz
2.0M -rw-r--r-- 1 stu411 stu411 2.0M 6月  14 00:29 time1_rep1_R1.fastq.gz
2.1M -rw-r--r-- 1 stu411 stu411 2.1M 6月  14 00:29 time1_rep1_R2.fastq.gz
2.0M -rw-r--r-- 1 stu411 stu411 2.0M 6月  14 00:29 time1_rep2_R1.fastq.gz
2.2M -rw-r--r-- 1 stu411 stu411 2.2M 6月  14 00:29 time1_rep2_R2.fastq.gz

#index文件
准备小鼠基因组的index文件：
下载地址：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz，下载后解压
最终路径：/lustre/home/acct-lctest/stu411/RNA-seq/index/mm10/genome

#注释文件
准备小鼠基因组的注释文件：
下载地址：ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/gencode.vM25.annotation.gtf.gz,下载后解压
最终路径：/lustre/home/acct-lctest/stu411/RNA-seq/annotation/gencode.vM25.annotation.gtf

########################################开始分析数据#######################################
学生从下面开始练习

########################################解压##############################################
gunzip fastqc time0_rep1_R1.fastq.gz
gunzip fastqc time0_rep1_R2.fastq.gz

#########################################质量控制########################################
fastqc time0_rep1_R1.fastq
fastqc time0_rep1_R2.fastq

#####################################reads比对#############################################
hisat2 -x /lustre/home/acct-lctest/stu411/RNA-seq/index/mm10/genome -1 time0_rep1_R1.fastq -2 time0_rep1_R2.fastq -S time0_rep1_align.sam 
samtools view -bS time0_rep1_align.sam > time0_rep1_unsorted.bam 
samtools sort time0_rep1_unsorted.bam -o time0_rep1_sorted.bam 

###等所有样本都跑完上述命令，运行下面命令

#####################################reads计数#############################################
featureCounts -p -t exon -g gene_id -a ./annotation/gencode.vM25.annotation.gtf -o sample_featurecount_txt ./time0_rep1_sorted.bam ./time0_rep2_sorted.bam ./time1_rep1_sorted.bam ./time1_rep2_sorted.bam 

#####################################准备DESeq2的输入文件##################################
#删除第一行
sed -i '1d' sample_featurecount_txt
#提取count
cut -f 1,7,8,9,10 sample_featurecount_txt > sample_featurecount_input.txt 


####################################差异表达分析###########################################
#进入R环境
module purge
module load miniconda3/4.7.12.1-gcc-4.8.5
source activate R
R

#DESeq2
BiocManager::install("DESeq2")
setwd('/lustre/home/acct-lctest/stu411/RNA-seq/process_data')
library("DESeq2")
input_data <- read.table("sample_featurecount_input.txt",header = TRUE,row.names =1)
input_data <- as.matrix(input_data)
condition <- factor(c(rep("time0",2),rep("time1",2)))
coldata <- data.frame(row.names=colnames(input_data),condition)
dds <- DESeqDataSetFromMatrix(countData =input_data,colData=coldata,design= ~condition)
#删除count数小于1或等于1的
dds1 <- dds[ rowSums(counts(dds)) > 1, ]
#进行差异基因分析
dds2 <- DESeq(dds1)
res_0.05 <- results(dds2,alpha=0.05)
res_0.05 <- res_0.05[order(res_0.05$padj),]
resdata_0.05 <- merge(as.data.frame(res_0.05),as.data.frame(counts(dds2,normalized=TRUE)),by="row.names",sort=FALSE)
write.csv(resdata_0.05,file="diff_gene_whole.csv",quote=F,row.names=F)
#按照自定义的阈值提取差异基因并导出
diff_gene_deseq2 <-subset(resdata_0.05, padj < 0.05 & abs(log2FoldChange) > 1)
write.csv(diff_gene_deseq2,file= "DEG.csv",quote=F,row.names=F)

#####################################功能富集分析#########################################
BiocManager::install("clusterProfiler")
library("clusterProfiler")
library("org.Mm.eg.db")
gene <- read.csv("DEG.csv",header = TRUE,sep=',')
ENSEMBL <- gsub("\\.\\d*", "", gene$Row.names) 
gene$Row.names <- ENSEMBL
#convert id
id <- gene[,1]
li <- bitr(id,fromType = "ENSEMBL",toType = c("GENENAME","SYMBOL","ENTREZID"),OrgDb =org.Mm.eg.db)
#go enrichment 
ego <- enrichGO(gene = li[,4],keyType = "ENTREZID",OrgDb = org.Mm.eg.db,pvalueCutoff = 0.05, pAdjustMethod = "BH",qvalueCutoff = 0.05,ont = "BP",readable = TRUE)
write.csv(as.data.frame(ego),"go_bp.csv",row.names =F)
#kegg enrichment 
ekegg <- enrichKEGG(li[,4], keyType = "kegg",organism = "mmu", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.05)
ekegg1 <- setReadable(ekegg, OrgDb = org.Mm.eg.db, keyType="ENTREZID")
write.csv(as.data.frame(ekegg1),"kegg.csv",row.names =F)

#########################################end##############################################


