STAR --runThreadN 6 --runMode genomeGenerate \
--genomeDir /home/sll/genome-sheep/Oar_rambouillet_v1.0-STAR \
--genomeFastaFiles GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
--sjdbGTFfile GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf \
--sjdbOverhang 149
--runThreadN： 线程数。
--runMode genomeGenerate： 构建基因组索引。
--genomeDir： 索引目录。（需提前建好）
--genomeFastaFiles： 基因组文件。
--sjdbGTFfile： 基因组注释文件。
--sjdbOverhang： reads长度减1
STAR --runThreadN 10 --genomeDir /home/sll/genome-sheep/Oar_rambouillet_v1.0-STAR \
--readFilesIn SRR17709920_1.filter.fastq.gz SRR17709920_2.filter.fastq.gz \
--readFilesCommand gunzip -c \
--outFileNamePrefix ./starmap.bam/SRR17709920
--readFilesIn ： paired reads文件
--outSAMtype ： 表示输出默认排序的bam文件，类似于samtools sort（还有--outSAMtype BAM Unsorted和--outSAMtype BAM Unsorted SortedByCoordinate）
--outFileNamePrefix ： 输出文件路径即前缀
- 常规比对
STAR --runThreadN 10 --genomeDir /home/sll/genome-sheep/Oar_rambouillet_v1.0-STAR \
--readFilesIn SRR17709920_1.filter.fastq.gz SRR17709920_2.filter.fastq.gz \
--readFilesCommand gunzip -c \
--outFileNamePrefix ./starmap.bam/SRR17709920
--readFilesIn ： paired reads文件
--outSAMtype ： 表示输出默认排序的bam文件，类似于samtools sort（还有--outSAMtype BAM Unsorted和--outSAMtype BAM Unsorted SortedByCoordinate）
--outFileNamePrefix ： 输出文件路径即前缀
- 建立样品索引，--sjdbFileChrStartEnd参数将所有样品的SJ.out.tab文件作为输入的annotated junction进行第二次建index。
STAR --runThreadN 20 --runMode genomeGenerate \
--genomeDir ./index \
--genomeFastaFiles /home/sll/genome-sheep/Oar_rambouillet_v1.0-STAR/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
--sjdbGTFfile /home/sll/genome-sheep/Oar_rambouillet_v1.0-STAR/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf \
--sjdbFileChrStartEnd *.out.tab --sjdbOverhang 149
./index 新建的索引目录
- 第二次比对，用第二次建立的index再一次对每个样品进行STAR比对
STAR --runThreadN 10 --genomeDir ./starmap.bam/index \
--readFilesIn SRR17709920_1.filter.fastq.gz SRR17709920_2.filter.fastq.gz \
--readFilesCommand gunzip -c \
--outFileNamePrefix ./starmap.bam/2_pass/SRR17709920_2pass
--readFilesCommand gunzip -c 如果使用的是.gz文件，则加此参数
######注：常规的差异分析、可变剪接分析建议使用HISAT2软件，这破软件好烦，如果是要用RNA数据call SNP那就使用它的2-pass模式，毕竟是GATK官网推荐的软件。
单样品示例
java -jar -Xmx15g /home/software/picard/build/libs/picard.jar AddOrReplaceReadGroups I=SRR17709911_2passAligned.out.sam \
O=./sample_RNA_SNP_calling/SRR17709911_sort_added.bam \
SO=coordinate \
RGID=SRR17709911 \
RGLB=mRNA \
RGPL=illumina \
RGPU=NovaSeq6000 \
RGSM=SRR17709911
全部样品一起（PS:这一步可以一起跑，不太占内存,）
ls *sort.bam | cut -d"_" -f 1 | sort -u | while read id ;
do (nohup java -jar -Xmx15g /home/software/picard/build/libs/picard.jar AddOrReplaceReadGroups I=${id}_sort.bam \
O=./sample_RNA_SNP_calling/${id}_sort_added.bam \
SO=coordinate \
RGID=${id} \
RGLB=mRNA \
RGPL=illumina \
RGPU=NovaSeq6000 \
RGSM=${id} &) ;
done
单样品示例
java -jar /home/software/picard/build/libs/picard.jar MarkDuplicates I=SRR17709911_sort_added.bam \
O=SRR17709911.sorted.addhead.markdup.bam \
M=SRR17709911.sorted.addhead.markdup.metrics.txt \
CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT
全部样品一起（PS:悠着点，一般64G运存的一次跑4个就行了，不然服务器崩了别怪我）
ls *sort_added.bam | cut -d"_" -f 1 | sort -u | while read id ;
do (nohup java -jar /home/software/picard/build/libs/picard.jar MarkDuplicates I=${id}_sort_added.bam \
O=${id}.sorted.addhead.markdup.bam M=${id}.sorted.addhead.markdup.metrics.txt \
CREATE_INDEX=true \
VALIDATION_STRINGENCY=SILENT &) ;
done
-- Split N trim and reassign mapping qualities不是单纯提取splitbam而是将所有bam转化成GATK可识别的格式- 这一步是
RNA-seq特异性的一步。因为mRNA转录本是主要由DNA的外显子exon可变剪切组合而成
单样品示例
/home/software/gatk-4.1.4.0/gatk SplitNCigarReads -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-I SRR17709911.sorted.addhead.markdup.bam \
-O SRR17709911.sorted.addhead.markdup.split.bam
全部样品一起（PS:一次最多4个）
ls *markdup.bam | cut -d"." -f 1 | sort -u | while read id ;
do (nohup /home/software/gatk-4.1.4.0/gatk SplitNCigarReads -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-I ${id}.sorted.addhead.markdup.bam \
-O ${id}.sorted.addhead.markdup.split.bam &) ;
done
单样品示例
/home/software/gatk-4.1.4.0/gatk HaplotypeCaller -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
--output-mode EMIT_ALL_CONFIDENT_SITES \
--dont-use-soft-clipped-bases \
-stand-call-conf 20 \
-ERC GVCF \
--max-mnp-distance 0 \
-I SRR17709911.sorted.addhead.markdup.split.bam \
-O SRR17709911.gvcf.gz
--output-mode EMIT_ALL_CONFIDENT_SITES
--dont-use-soft-clipped-bases 表示GATK不考虑比对时产生的soft clipped的碱基，以减少假阳性和假阴性的变异
-stand-call-conf 相当于一个可信度打分，转录的默认是20，全基因组会考虑用30
-ERC GVCF 输出gvcf方便后续的合并
--max-mnp-distance 0为不输出mnp,否则后续合并时会报错，但是一般默认为0，就挺奇怪的
全部样品一起（PS:一次最多4个）
ls *markdup.split.bam | cut -d"." -f 1 | sort -u | while read id ;
do (nohup /home/software/gatk-4.1.4.0/gatk HaplotypeCaller -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
--output-mode EMIT_ALL_CONFIDENT_SITES \
-ERC GVCF \
--dont-use-soft-clipped-bases \
-stand-call-conf 20 \
-I ${id}.sorted.addhead.markdup.split.bam \
-O ${id}.gvcf.gz &) ;
done
！！！如果在合并过程中出现报错
“A USER ERROR has occurred: Bad input: Combining gVCFs containing MNPs is not supported. Unknown contained a MNP at NC_037328.1:994538”
则可以在HaplotypeCaller时添加参数--max-mnp-distance 0
nohup /home/software/gatk-4.1.4.0/gatk CombineGVCFs -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-V SRR17709910.gvcf.gz \
-V SRR17709911.gvcf.gz \
-V SRR17709912.gvcf.gz \
-O combined.vcf &
-V 要合并的gvcf文件
-O 输出的vcf文件名称
/home/software/gatk-4.1.4.0/gatk GenotypeGVCFs -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-V combined.vcf \
-O combined.genotype.vcf
/home/software/gatk-4.1.4.0/gatk SelectVariants --select-type SNP \
-V combined.genotype.vcf \
-O combined.genotype.snp.vcf
RNA数据不支持VQSR,使用 FS > 30.0 & QD < 2.0，GATK建议过滤掉在35bp范围内出现3个以上的SNP的情况（-window 35 -cluster 3）推荐命令
/home/software/gatk-4.1.4.0/gatk VariantFiltration -R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-V combined.genotype.snp.vcf \
-O combined.genotype.snp.filter.vcf \
--window 35 --cluster 3 \
--filter-name "FS>30.0" --filter-expression "FS>30.0" \
--filter-name "QD<2.0" --filter-expression "QD<2.0"
/home/software/gatk-4.1.4.0/gatk SelectVariants --exclude-filtered \
-R /home/sheep-reference/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.fna \
-V combined.genotype.snp.filter.vcf \
-O combined.genotype.snp.filter.retain.vcf
# 单个文件
wget -c -t 0 -O SRR13107018.sra https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR13107018/SRR13107018
#-c -t 配合使用可以防止下载数据的过程中链接中断的问题，-O则可以指定下载路径和文件名。
# 多个文件一起下
/home/software/sratoolkit.2.9.6-1-centos_linux64/bin/prefetch -O ./ --option-file SRR_Acc_List.txt
ls *sra | while read id ;
do (nohup /home/software/sratoolkit.2.9.6-1-centos_linux64/bin/fastq-dump --gzip --split-3 $id &) ;
done
ls *fastq.gz | cut -d"_" -f 1 |sort -u | while read id ;
do (nohup fastp -i ${id}_1.fastq.gz -I ${id}_2.fastq.gz -g -q 5 -n 15 -l 150 -u 50 -o ${id}_1.filter.fastq.gz -O ${id}_2.filter.fastq.gz -h ${id}.fastp.html &) ;
done
# 去除带接头(adapter)的双端读长(默认开启，可用-A关闭)；
# 去掉单端读长中含 N(N 表示无法确定碱基信息)的碱基比例大于10%部分（默认开启）
# -q 设置碱基质量阈值，默认是15，phred质量评分≥Q15
# -n 一条read中N碱基的数量超过了多少个，那么这条read就被删除，默认是5（即这条read里有5个N）
# -g 进行PolyG尾的去除
# -u 允许百分之多少的碱基不合格（0-100），默认是40（就是说40%），超过这个比例，整条read就被删除
# -l read小于这个参数设定长度会被丢弃或删除，默认是15,由你的测序bp决定
# -j, --json 输出的json报告文件名，以“.json”结尾
# -h, --html 输出的html报告文件名，以“.html”结尾
由于使用hisat2进行序列比对，因此也用它建索引
## 下载基因组序列
mkdir -p genome/ARS-UCD1.2 && cd genome/ARS-UCD1.2
nohup wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/263/795/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz &
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz
## 下载gtf文件并解压
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Bos_taurus/latest_assembly_versions/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
## 建索引（hisat2软件）
mkdir index & cd index
nohup hisat2-build -p 4 GCF_002263795.1_ARS-UCD1.2_genomic.fna GCF_002263795.1_ARS-UCD1.2_genomic > hisat2.log &
mkdir hismap.sam
ls *filter.fastq.gz|cut -d"_" -f 1 |sort -u | while read id ;
do (nohup hisat2 -p 8 -x /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic -1 ${id}_1.filter.fastq.gz -2 ${id}_2.filter.fastq.gz -S hismap.sam/${id}.hismap.sam &) ;
done
#绵羊参考基因组
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic
#牛参考基因组
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic
# -p: 线程数
# -x: 基因组索引前缀。所下的基因组索引为多个文件，索引前缀到genome为止。
# -1/-2: fastq输入文件。当输入为单端测序时使用-U 指定输入。单端或双端的输入都可使用逗号分隔输入多个样本，或使用多次-1 -2 / -U 指定多个输入。e.g.: -U sample1.fq.gz,sample2.fq.gz -U sample3.fq.gz
# -S: 输出sam文件路径。
cd hisout.sam
ls *sam | cut -d"." -f 1 |sort -u | while read id ;
do (nohup samtools view -S ${id}.hismap.sam -b > ${id}.hismap.bam &) ;
done
ls *hismap.bam | cut -d"." -f 1 |sort -u | while read id ;
do (nohup samtools sort -@ 8 ${id}.hismap.bam -o ${id}_sort.bam &) ;
done
# sort: 进行排序
# -@: 线程数
# -o: 输出bam文件
# 最后一项为输入文件
ls *sort.bam | cut -d"_" -f 1 | sort -u | while read id ;
do (nohup samtools index ${id}_sort.bam ${id}_sort.bam.index &) ;
done
## 使用samtools index建立索引，使之快速访问bam文件
nohup featureCounts -p -t exon -g gene_id -M -T 8 \
-a /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf \
-o all.featurecounts.txt \
*_sort.bam &
#绵羊gtf文件
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf
#牛gtf文件
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf
# featureCounts 进行定量，统计比对在这个基因的坐标上的read数
# -t 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”
# -p 双端数据时需要
# -a 转录组注释文件
# -o 输出文件
##报错处理
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id | wc -w #检查原始gtf注释文件gene_id个数
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w #检查空gene_id值的行数
sed -i -e '/gene_id\ \"\"\;/d' GCF_002263795.1_ARS-UCD1.2_genomic.gtf #删除包括空gene_id值的行
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w #修改后的检查一下 检查包括空gene_id值的行 0为已删除掉
multiqc all.featurecounts.txt.summary -o all.counts.summary
awk -F '\t' '{print $1,$7,$8,$9,$10,$11,$12}' OFS='\t' all.featurecounts.txt > all_fcount.matrix.txt
# 1列为基因id，7列及以后为样本定量信息
# \t 表示以制表符分割开来
# 安装包
#####################################################
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
install.packages("corrplot")
install.packages("PerformanceAnalytics")
install.packages("factoextra")
#####################################################
# 1 数据读入和处理
rm(list=ls())
countdata<- read.table("all_fcount.matrix.txt", header=TRUE,row.names = 1) # 导入数据
head(countdata) # 查看数据前几行(列名 太长自己展示看)
countdata.filter<-countdata[rowSums(countdata)>=0&apply(countdata,1,function(x){all(x>=10)}),] # 过滤低丰度基因，过滤featurecounts后，每个样本计数小于等于10的！！！！！！！
head(countdata.filter)
colnames(countdata.filter) <- c("RIBEN_1","RIBEN_2","RIBEN_3","CHINA_1","CHINA_2","CHINA_3")# 修改列名
head(countdata.filter)
dim(countdata.filter)
#####################################################
# 2 构建dds矩阵
library(DESeq2)
condition <- factor(c("RIBEN","RIBEN","RIBEN","CHINA","CHINA","CHINA")) # 赋值因子，即变量
coldata <- data.frame(row.names=colnames(countdata.filter), condition) # 创建一个condition数据框
dds <- DESeqDataSetFromMatrix(countData=countdata.filter, colData=coldata, design=~condition)# 构建dds矩阵(后面很多分析都是基于这个dds矩阵)
#####################################################
# 3 DESeq2进行差异分析
dds <- DESeq(dds) # DESeq分析
resdata<- results(dds,contrast = c("condition","RIBEN","CHINA")) # 此为"RIBEN"比"CHINA"
table(resdata$padj<0.05 ) # Benjamini-Hochberg矫正后的p<0.05的基因数！！！！！！！
res_padj <- resdata[order(resdata$padj), ] # 按照padj(矫正后的p值)列的值排序
write.table(res_padj,"diffexpr_padj_results.txt",quote = F,sep = '\t') # 将结果文件保存到本地，打开在第一列头加gene
####################################################
# 4 DESeq2分析中得到的resdata进行绘制火山图
rm(list=ls())
resdata <- read.table("diffexpr_padj_results.txt",header = T,sep = '\t',row.names = 1) # 加载DESeq2中生成的resdata文件
resdata$label <- c(rownames(resdata)[1:10] ,rep(NA,(nrow(resdata)-10))) # 对前10个基因进行标注，看自己需求，也可不要这一步
library(ggplot2)
ggplot(resdata,aes(x=log2FoldChange,y=-log10(padj))) +
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "#999999")+ # 横向水平参考线，显著性 --p值
geom_vline(xintercept = c(-0.585 , 0.585),linetype = "dashed", color = "#999999")+ # 纵向垂直参考线，差异倍数 --foldchange，此处为1.5倍
geom_point(aes(size = -log10(padj),color = -log10(padj))) + # 散点图
scale_color_gradientn(values = seq(0,1,0.2),
colors = c("#39489f","#39bbec","#f9ed36","#f38466","#b81f25"))+ # 指定颜色渐变模式
scale_size_continuous(range = c(0,3))+ # 指定散点渐变模式
## 主题调整
theme_bw()+ # theme_grey()默认主题，theme_bw()白色背景，theme_classic()经典主题
theme(panel.grid.major = element_blank(), # 删除主网格线
panel.grid.minor = element_blank(), # 删除次网格线
legend.position = c(0.08,0.9), # 设置图例位置
legend.justification = c(0,1))+
guides(col = guide_colourbar(title = "-Log10_q-value"),
size = "none")+ # 设置部分图例不显示，size = "none"
geom_text(aes(label=c(label),color = -log10(padj)), size = 3, vjust = 1.5, hjust = 1)+ # 添加标签,也就是上述resdata$label，如不需要，可删除此列代码
xlab("Log2FC")+ylab("-Log10(FDR q-value)") # 修改坐标轴名称
ggsave("vocanol_plot.pdf", height = 9 , width = 10)
####################################################
# 5 筛选上下调差异基因
subset(resdata,pvalue < 0.05) -> diff # 先筛选pvalue < 0.05的行！！！！！
subset(diff,log2FoldChange < -0.585) -> down # 筛选出下调的,此处定位1.5倍
subset(diff,log2FoldChange > 0.58) -> up # 筛选出上调的
dim(down) # 查看数据维度，即几行几列
dim(up)
up_names<-rownames(up) # 导出上、下调基因的那一列
write.table(up_names,'up_gene.txt',quote = F,sep = '\t',row.names = F) # 将结果写到本地
down_names <- rownames(down)
write.table(down_names,'down_gene.txt',quote = F,sep = '\t',row.names = F)
####################################################
# 6 表达矩阵探索,选取差异表达的基因做热图 (注意是用矫正后p值排完序的)
library(pheatmap)
choose_gene=head(rownames(resdata),50) # 定义前50个为筛选的基因
choose_matrix=countdata.filter[choose_gene,] # 抽取差异表达显著的前50个基因
choose_matrix=t(scale(t(choose_matrix))) # 用t函数转置，scale函数标准化
pheatmap(choose_matrix)
####################################################
# 7 PCA分析 使用对dds矩阵处理后的vst或rld值
dds_pca <- estimateSizeFactors(dds) # 计算每个样本的归一化系数
vsd <- vst(dds_pca,blind = FALSE) # 方差稳定变换
rld <- rlog(dds_pca,blind = FALSE) # 正则化对数变换
# rlog函数将count data转换为log2尺度，以最小化有small counts的行的样本间差异，并使library size标准化
# vst函数快速估计离散趋势并应用方差稳定变换。
# 数据集小于30个样品可以用rlog，数据集大于30个样品用vst，因为rlog速度慢
# 转换的目的是，为了确保所有基因有大致相同的贡献
plotPCA(vsd, intgroup=c("condition"))
plotPCA(rld, intgroup=c("condition"))
# 安装包
#####################################################
install.packages("BiocManager")
BiocManager::install(version = "3.13")
BiocManager::install("gprofiler2")
BiocManager::install("clusterProfiler")
BiocManager::install("AnnotationHub")
#####################################################
# 注释库的安装
## 方法1 在线
library(AnnotationDbi)
library(AnnotationHub)
hub = AnnotationHub()
hub$species[which(hub$species == "Ovis aries")] # 寻找AnnotationHub中是否有目标物种的注释库
org <- hub[hub$rdataclass == "OrgDb",] # 指定OrgDb数据库格式，而不是其他数据库格式
query(org, "Ovis aries") # 查找目标物种的数据库，找.sqlite文件对应的AH号
sheep.db <- hub[['AH6429']] # 下载注释库为sheep.db，注释库加载完成
columns(sheep.db)
## 方法2 本地
setwd("D:/生信/sheep.db") # 下载的数据库保存目录，加载对应目录
sheep.db <- AnnotationDbi::loadDb('org.Ovis_aries.eg.sqlite') # 在Bioconder下载物种对应的注释库.sqlite，加载完成
columns(sheep.db)
## 牛
setwd("D:/生信/cattle.db")
bos.db <- AnnotationDbi::loadDb('org.Bt.eg.db.sqlite')
columns(bos.db)
#####################################################
# 1 GO分析（上下调基因分开做，可用于BP,CC,MF分开画图）
library(clusterProfiler)
genes <- read.table('GENEID.txt', header = T, stringsAsFactors = FALSE,sep = '\t')# 读取基因列表文件中的基因名称，GENEID.txt为一行的无表头的基因ID
enrich.go <- enrichGO(gene = genes$V1, # 基因列表文件中的基因名称
OrgDb = bos.db, # 指定物种的基因数据库
keyType = 'SYMBOL', # 指定给定的基因名称类型，例如这里以 SYMBOL 为例
ont = 'CC', # 可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
pAdjustMethod = 'BH', # 指定 p 值校正方法
pvalueCutoff = 1, # 指定 p 值阈值，不显著的值将不显示在结果中
qvalueCutoff = 1 # 指定 q 值阈值，不显著的值将不显示在结果中
)
# 例如上述指定 ALL 同时计算 BP、MF、CC，这里将它们作个拆分后输出
BP <- enrich.go[enrich.go$ONTOLOGY=='BP', ]
CC <- enrich.go[enrich.go$ONTOLOGY=='CC', ]
MF <- enrich.go[enrich.go$ONTOLOGY=='MF', ]
write.table(as.data.frame(BP), 'HUgo.BP.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(CC), 'HUgo.CC.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(MF), 'HUgo.MF.txt', sep = '\t', row.names = FALSE, quote = FALSE)
# 2 GO分析（上下调一起做，可以看整体结果，并可选取其中的部分画图）
gene_down <- read.table('up_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
gene_up <- read.table('down_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
library(clusterProfiler)
deg.up.gene <- bitr(gene_up$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID # 转基因ID为ENTREZID
deg.down.gene <- bitr(gene_down$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
deg.ei.ls <- list(up = deg.up.gene, down = deg.down.gene) # 将up 和 down基因合并到一个表格
library("org.Bt.eg.db")
GO.cmp <- compareCluster(deg.ei.ls,
fun = "enrichGO",
qvalueCutoff = 0.1,
pvalueCutoff = 0.01,
pAdjustMethod = 'BH',
ont = 'BP', # 可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者，建议分开做输出，默认按照pvalue升序排列的
OrgDb = bos.db)
write.table(GO.cmp@compareClusterResult,'GO-up-down.txt',sep = '\t',quote = F) # 结果写到本地GO-up-down.txt
#####################################################
# 3 GO 条形图（上下调一起，不过BP,CC,MF等需要分开做）
library(ggplot2)
library(ggsci)
ddf<-read.table('GO-up-down.txt',header = T,sep = '\t') # 上一步输出的上下调富集结果，一起做的那个
dat=ddf
dat$Cluster = ifelse(dat$Cluster =='up',1,-1)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue <- dat$pvalue*dat$Cluster
dat=dat[order(dat$pvalue,decreasing = F),]
ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)),
y=pvalue,fill=Cluster)) +
geom_bar(stat="identity",width = 0.78,position = position_dodge(0.7)) +
scale_fill_gradient(low="#3685af",high="#af4236",guide = FALSE)+
scale_x_discrete(name ="GO Terms") +
#scale_y_discrete(labels =Description )
scale_y_continuous(name ="-log10Pvalue",expand = c(0,0)) +
coord_flip() + # 翻转坐标
theme(panel.grid.major.y = element_blank(),panel.grid.minor.y = element_blank())+ # 删除纵向网络
theme(axis.text.y = element_blank())+ # 刻度标签空白
theme(axis.ticks = element_blank())+ # 刻度为空白
theme(panel.border = element_blank())+ # 边框为空白
theme(axis.text=element_text(face = "bold",size = 15),
axis.title = element_text(face = 'bold',size = 15), # 坐标轴字体
plot.title = element_text(size = 20,hjust = 0.3),
legend.position = "top",panel.grid = element_line(colour = 'white'))+
geom_text(aes(label=Description),size=3.5,hjust=ifelse(dat$pvalue>0,1.5,-0.5))+ # 两侧加上标签文字
ggtitle("The most enriched GO-BP terms") # 添加标题
# 4 GO气泡图（可上下调一起做，也不用将BP,CC,MF分开，最好选择适当数量，KEGG同样适用）
BiocManager::install("openxlsx")
library(ggplot2)
library(openxlsx) # 这个包是用来导入格式为XLSX的数据的
goinput <- read.xlsx("go-barplot.xlsx") # 载入GO富集分析的结果文件
# 设置XY轴
x=goinput$GeneRatio # 设置以哪个数据为X轴
y=factor(goinput$Description,levels = goinput$Description) # 设置Y轴
p = ggplot(goinput,aes(x,y)) # 开始以X轴和Y轴绘制画布
p
p1 = p + geom_point(aes(size=Count,color=-0.5*log(pvalue),shape=ONTOLOGY,))+
scale_color_gradient(low = "SpringGreen", high = "DeepPink") # 以count为点的大小，以颜色的由嫩绿到深粉红代表P值的对数值，以点的性状代表GO类型（MF,CC,BP）
p1
p2 = p1 + labs(color=expression(-log[10](pvalue)),
size="Count",
x="Ratio", # 设置x轴标签
y="Go_term", # 设置x轴标签
title="Go enrichment of test Genes") # 设置图标题
p2
#####################################################
# 5 KEGG, 上下调一起做（适合看上下调整体的结果）
gene_down <- read.table('up_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
gene_up <- read.table('down_gene.txt', header = T, stringsAsFactors = FALSE,sep = '\t')
library(clusterProfiler)
deg.up.gene <- bitr(gene_up$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
deg.down.gene <- bitr(gene_down$x, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = bos.db)$ENTREZID
deg.ei.ls <- list(up = deg.up.gene, down = deg.down.gene)
search_kegg_organism("taurus", by="scientific_name") # 寻找对应动物的organism
library(clusterProfiler)
R.utils::setOption("clusterProfiler.download.method",'auto') # 一定要用这个，不然报错别怪我哦
kegg.cmp <- compareCluster(deg.ei.ls,
fun = "enrichKEGG",
qvalueCutoff = 1,
pvalueCutoff = 1,
pAdjustMethod = 'BH',
organism = "bta")
write.table(kegg.cmp@compareClusterResult,'kegg-up-down.txt',sep = '\t',quote = F)
# 6 kegg， 上下调分开做（适合画图）
library(clusterProfiler)
genes <- read.table('txt', header = T, stringsAsFactors = FALSE,sep = '\t')
covID <- bitr(genes$V1, fromType = "SYMBOL",
toType="ENTREZID",OrgDb= bos.db , drop = TRUE) # 转换为entrezID，KEGG只识别这种ID，报错请检查自己输入的文件，对自己的代码有自信
write.table(covID,'ENTREZID-up.txt',sep = '\t',quote = F) # 写出去保存,也可以不保存，看个人需求
#每次打开R计算时，它会自动连接kegg官网获得最近的物种注释信息，因此数据库一定都是最新的
search_kegg_organism("taurus", by="scientific_name") # 寻找对应动物的organism
R.utils::setOption("clusterProfiler.download.method",'auto')
enrich_kegg <- enrichKEGG(gene = covID$ENTREZID, # 基因列表文件中的基因名
organism = 'bta', # 指定物种，牛为bta
keyType = 'kegg', # kegg富集
pAdjustMethod = 'BH', # 指定矫正方法
pvalueCutoff = 1, # 指定p值阈值，不显著的值将不显示在结果中
qvalueCutoff = 1) # 指定q值阈值，不显著的值将不显示在结果中
write.table(enrich_kegg@result,'kegg-down.txt',quote = F,sep = '\t')
#####################################################
# 7 KEGG富集气泡图(上下调分开做推荐)
library(ggplot2)
dotplot(enrich_kegg,
showCategory = 15, # 展示前15个点，默认为10个
color = 'pvalue',
title = "kegg_dotplot")
# 8 KEGG富集条形图
barplot(enrich_kegg,
showCategory = 15,
color = 'pvalue',
title = "kegg_dotplot")
#####################################################
# GO气泡图
# BiocManager::install("openxlsx")
library(ggplot2)
library(openxlsx)
goinput <- read.xlsx("go.xlsx") # 载入GO富集分析的结果文件
x=goinput$GeneRatio # 设置以哪个数据为X轴
y=factor(goinput$Description,levels = goinput$Description) # 设置Y轴
p = ggplot(goinput,aes(x,y)) # 开始以X轴和Y轴绘制画布
p1 = p + geom_point(aes(size=Count,color=-0.5*log(pvalue),shape=ONTOLOGY,))+
scale_color_gradient(low = "SpringGreen", high = "OrangeRed")
p1
p1 + labs(color=expression(-log[10](pvalue)),
size="Count",
x="GeneRatio",
y="Go_term",
title="Go enrichment of test Genes")
# KEGG气泡图
library(ggplot2)
library(openxlsx)
kegginput <- read.xlsx("kegg.xlsx")
x=kegginput$pvalue
y=factor(kegginput$Term,levels = kegginput$Term)
p = ggplot(kegginput,aes(x,y))
p1 = p + geom_point(aes(size=Count,color=-0.5*log(pvalue)))+
scale_color_gradient(low = "BLUE", high = "OrangeRed") # 以颜色的由暗蓝到橙红代表P值的对数值
p2 = p1 + labs(color=expression(-log[10](pvalue)),
size="Count",
x="P value",
y="KEGG",
title="KEGG pathway of Target Genes")
p2
# FPKM的计算
# 1 数据载入
rm(list=ls())
options(stringsAsFactors = F)
library(tidyverse)
library(data.table) # 可多核读取文件
a1 <- fread('all.featurecounts.txt', header = T, data.table = F) # 载入featurecounts后产生的row count文件，第一列设置为列名
# 2 counts矩阵的构建
counts <- a1[,7:ncol(a1)] # 截取样本基因表达量的counts部分作为counts
rownames(counts) <- a1$Geneid # 将基因名作为行名
geneid_efflen <- subset(a1,select = c("Geneid","Length")) # 从featurecounts 原始输出文件counts.txt中提取Geneid、Length(转录本长度)
colnames(geneid_efflen) <- c("geneid","efflen")
geneid_efflen_fc <- geneid_efflen
dim(geneid_efflen)
efflen <- geneid_efflen[match(rownames(counts),
geneid_efflen$geneid),"efflen"] # 取出counts中geneid的对应的efflen
# 3 FPKM/RPKM (Fragments/Reads Per Kilobase Million ) 每千个碱基的转录每百万映射读取的Fragments/reads
# RPKM与FPKM分别针对单端与双端测序而言，计算公式是一样的
counts2FPKM <- function(count=count, efflength=efflen){
PMSC_counts <- sum(count)/1e6 # counts的每百万缩放因子(“per million” scaling factor)深度标准化
FPM <- count/PMSC_counts # 每百万reads/Fragments (Reads/Fragments Per Million) 长度标准化
FPM/(efflength/1000)
}
FPKM <- as.data.frame(apply(counts,2,counts2FPKM))
colnames(FPKM) <- c("Simmental_1","Simmental_2","Simmental_3","Wagyu_1","Wagyu_2","Wagyu_3") # 修改列名，也可不修改，按照自己需求
FPKM <- FPKM[rowSums(FPKM)>=1,] # 去除全部为0的列
colSums(FPKM)
# TPM 当前推荐使用 TPM 进行相关性分析、PCA分析等 (Transcripts Per Kilobase Million) 每千个碱基的转录每百万映射读取的Transcripts
counts2TPM <- function(count=count, efflength=efflen){
RPK <- count/(efflength/1000) # 每千碱基reads (reads per kilobase) 长度标准化
PMSC_rpk <- sum(RPK)/1e6 # RPK的每百万缩放因子(“per million” scaling factor )深度标准化
RPK/PMSC_rpk
}
TPM <- as.data.frame(apply(counts,2,counts2TPM))
colnames(TPM) <- c("Zebu_1","Zebu_2","Zebu_3","Zebu_4","Zebu_5","Holstein_1","Holstein_2","Holstein_3","Holstein_4","Holstein_5")
TPM <- TPM[rowSums(TPM)>0,]
colSums(TPM)
-
可变剪切（differential splicing）也叫做选择性剪切alternative splicing, 指的是在mRNA前体到成熟mRNA的过程当中，不同的剪切方式使得同一个基因可以产生多个不同的成熟mRNA, 最终产生不同的蛋白质 https://zhuanlan.zhihu.com/p/409865441
-
rMATs软件可以识别5种可变剪接事件：Skippedexon(SE) 外显子跳跃、Alternative5’ splice site(A5SS) 5’端可变剪切、Alternative3’ splice site(A3SS) 3’端可变剪切、Mutuallyexclusive exons(MXE) 互斥可变外显子、Retainedintron(RI) 内含子保留
可以用conda直接安装
conda install -c bioconda rmats
rmats.py --b1 b1.txt \
--b2 b2.txt \
--gtf /home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf \
--od AS \
--tmp tmp \
-t paired \
--readLength 150 \
--cstat 0.0001 \
--nthread 10
--b1 b1.txt 输入sample1的txt格式的文件，文件内以逗号分隔重复样本的bam文件名
--b2 b2.txt 输入sample2的txt格式的文件，文件内以逗号分隔重复样本的bam文件名
-t readType 双端测序则readType为paired，单端测序则为single
--readLength 测序reads的长度,可以从质控报告看
--gtf gtfFile 需要输入的gtf文件
--od outDir 所有输出文件的路径（文件夹）
--nthread 设置线程数
--cstat The cutoff splicing difference. The cutoff used in the null hypothesis test for differential splicing
--statoff 进行单样本或者是单组的分析，并跳过统计分析
rmats2sashimiplot --b1 SRR17709921_sort.bam,SRR17709920_sort.bam,SRR17709917_sort.bam \
--b2 SRR17709910_sort.bam,SRR17709918_sort.bam,SRR17709919_sort.bam \
-t SE \
-e SE.MATS.JC.txt \
--l1 DP_L \
--l2 Han_L \
-o SE_plot &
# 可以将需要可视化的基因进行筛选，重新做成SE.MATS.JC.txt这种文件，然后可视化就可以了
rmats2sashimiplot --b1 SRR17709911_sort.bam,SRR17709912_sort.bam,SRR17709913_sort.bam \
--b2 SRR17709916_sort.bam,SRR17709915_sort.bam,SRR17709914_sort.bam \
-t SE \
-e SE.MATS.JC.txt \
--l1 DP_M \
--l2 Han_M \
-o M_SE_plot
--b1 B1 sample_1 in bam format(s1_rep1.bam[,s1_rep2.bam])
--b2 B2 sample_2 in bam format(s2_rep1.bam[,s2_rep2.bam])
-t rMATS结果中产生的可变剪切类型{SE,A5SS,A3SS,MXE,RI}
-e EVENTS_FILE The rMATS output event file (Onlyif using rMATSformat result as event file).
--l1 L1 The label for first sample.
--l2 L2 The label for second sample.-o OUT_DIR The output directory.
######会输出好几种文件，其中.MATS.JC.txt是我们要用到的文件
ID GeneID geneSymbol chr strand 1stExonStart_0base 1stExonEnd 2ndExonStart_0base 2ndExonEnd upstreamES upstreamEE downstreamES downstreamEE ID IJC_SAMPLE_1 SJC_SAMPLE_1 IJC_SAMPLE_2 SJC_SAMPLE_2 IncFormLen SkipFormLen PValue FDR IncLevel1 IncLevel2 IncLevelDifference
0 "MS.gene23798" NA chr8.4 - 30758609 30758704 30759122 30759182 30758025 30758095 30760455 30760527 0 1 11 7 9 209 244 0.0120878457309 0.0604392286545 0.096 0.476 -0.38
1 "MS.gene61989" NA chr7.2 - 80697619 80697769 80704270 80704420 80697113 80697232 80705567 80706851 1 1 1 0 3 298 298 0.102057409464 0.34019136488 0.5 0.0 0.5
●ID: 官网描述“rMATS event id”，其实就是序号
●GenelD: 可变剪接事件所在基因编号
●geneSymbol: 可变剪接事件所在基因名称
●chr: 可变剪接事件所在染色体
●strand: 可变剪接事件所在染色体链的方向
●1stExonStart_0base: 第一个可变剪接事件跳跃外显子的起始位置，以0开始计数
●1stExonEnd: 第一个可变剪接事件跳跃外显子的终止位置
●2ndExonStart_0base:第二个可变剪接事件跳跃外显子的起始位置，以0开始计数
●2ndExonEnd: 第二个可变剪接事件跳跃外显子的终止位置
●upstreamES: 可变剪接事件跳跃外显子的上游exon起始位置
●upstreamEE: 可变剪接事件跳跃外显子的上游exon终止位置
●downstreamES: 可变剪接事件跳跃外显子的下游exon起始位置
●downstreamEE: 可变剪接事件跳跃外显子的下游exon终止位置
●ID: 同上
●IJC_SAMPLE_1: 样本一在inclusion junction（IJC）下的count数，重复样本的结果以逗号分隔
●SJC_SAMPLE_1: 样本一在skipping junction（SJC）下的count数，重复样本的结果以逗号分隔
●IJC_SAMPLE_2: 样本二在inclusion junction（IJC）下的count数，重复样本的结果以逗号分隔
●SJC_SAMPLE_2: 样本二在skipping junction（SJC）下的count数，重复样本的结果以逗号分隔
●IncFormLen: 可变剪接事件Exon Inclusion Isoform的有效长度
●SkipFormLen: 可变剪接事件Exon Skipping Isoform的有效长度
●*PValue*: 两组样本间可变剪接事件表达差异显著性p值
●*FDR*: 可变剪接事件表达差异显著性FDR值
●*IncLevel1*: 处理组可变剪接事件Exon Inclusion Isoform在两个Isoform总表达量的比值，也就是**PSI**
●*IncLevel2*: 对照组可变剪接事件Exon Inclusion Isoform在两个Isoform总表达量的比值，也就是**PSI**
●*IncLevelDifference*: IncLevel1与IncLevel2的差值，和**dPSI（different percent spliced in）**差不多
