鉴定与发酵分枝杆菌感染介导的坏死神经元细胞死亡相关的RNA-Seq分析
本文作者：2020级生物信息国重班何思成
RNA-Seq基础分析流程参考：http://ftp.genek.cn:8888/20220919_mjvlyt/%E8%BD%AC%E5%BD%95%E7%BB%84%E5%9F%BA%E7%A1%80%E5%88%86%E6%9E%901014.pdf
参考的文献资料：Transcriptome analysis of necrotic neuronal cell death mediated by M. fermentans infection
Accession：PRJNA900032 ID: 900032
Data Type：Raw sequence reads
Submission：Registration date: 10-Nov-2022 Seoul national university
一. 相关研究背景
关于发酵支原体感染介导的坏死神经元细胞死亡相关的差异表达基因，相关NCBI上显示的文章并没有详细介绍，但研究人员经分离RNA构建文库，按照illumina测序后获得了三组测序库，其测序后，我将通过生物信息学方法和技术对三组RNA-Seq进行基因表达差异分析（三组技术重复，没有生物学重复）
这一部分是Linux的关键代码操作
SRA数据库显示的数据如下
image.png
下载原始数据sra
kingfisher get -r SRR22252649 -m ena-ftp aws-cp aws-http prefetch
下载参考基因组数据
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.gff.gz
处理SRA文件并解压
fastq-dump --split-files --gzip SRR22252649.sra
##因为是双端测序所以加上--split-files，gzip压缩处理更方便
质控过滤
fastp -i SRR222526${k}_1.fastq.gz -I SRR222526${k}_2.fastq.gz -o ../after-fastp/report${k}_1.fastq.gz -O ../after-fastp/report${k}_2.fastq.gz -h ../after-fastp/report${k}.html -j ../after-fastp/report${k}.json
处理参考基因组
seqtk seq -l 60 GCF_000001405.40_GRCh38.p14_genomic.fna.gz>gene.fa
##将参考基因组处理为每行60个碱基
gunzip GCF_000001405.40_GRCh38.p14_genomic.gff.gz
gffread -T -o genes.gtf GCF_000001405.40_GRCh38.p14_genomic.gff
##gff转为gtf
为参考基因组构建索引
nohup bwa index ../21.mapping/Gene.fa -p bwa_index &
##后台运行
比对到参考基因组,sam转为bam
samtools faidx Gene.fa
samtools view -bhS -t Gene.fa.fai -o Gene.bam BM1.sam
表达定量
nohup Rscript ~/hesic/10.software/RunFeatureCounts/run-featurecounts.R -b Gene.bam -g ~/hesic/rna/genes.gtf -f exon -a gene_id -o Gene_1 -s 2 &
#看一下log文件截图，count文件截图
合并成表达矩阵
ls ../count/use_count/*.count > genes.quantfiles.txt #将所有生成的*count文件路径打印到txt文件中，方便下一步操作
perl ../10.software/RunFeatureCounts/abundance_estimates_to_matrix.pl --est_method featureCounts --out_prefix gene --quant_files genes.quantfiles.txt
使用这个perl程序（这个是基因课老师教的，具体来源我也不知道）合并成表达矩阵
差异表达分析
perl /pub/software/trinityrnaseq-v2.11.0/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix ../gene.counts.matrix --method edgeR --dispersion 0.4 --samples_file ../hm_info.txt --contrasts ../contrast.txt
##没有生物学重复，自行设置dispersion 0.4
head -1 1v2gene.DE_results >./merge.DE_results
cat ./*.DE_results |grep -v "^sampleA" >>merge.DE_results
下面是Rstudio的关键代码操作
导入我的样本收集表
（其中hum_info.txt的内容是gene1 vs gene2,gene2 vs gene3, gene1 vs gene3）
library(tidyverse)
read.table(file = "~/hesic/hm_info.txt") %>%
rename( V1= "group", V2= "sample") -> sample_info
sample_info
导入表达矩阵
gene_exp = read.table(file = "~/hesic/gene.TMM.TPM.matrix")
简单查看相关系数（如下图）
sample_cor = round(cor(gene_exp),digits = 7)
sample_cor
image.png
绘制热图分析
library(pheatmap)
pheatmap(sample_cor,
cluster_cols = F,
cluster_rows = F,
cellwidth = 90,
cellheight = 90,
border_color= "white",
fontsize = 9,
angle_col = 90,
display_numbers = T,
fontsize_number = 4
)
image.png
查看基因距离层次
plot(hclust(dist(t(gene_exp))))
image.png
主成分分析（如下图）
library(PCAtools)
p1 = pca(gene_exp,removeVar = 0.3,metadata =sample_info)
screeplot(p1) ##查看各个成分之间的解释度
biplot(p1,
x="PC1",
y="PC2",
colby = "group",
legendPosition = "right",
hline = 1,vline = 1,
encircle = T)
image.png
差异基因探索
deg_files = list.files(path = "~/hesic/31_analysis/",
pattern = "DE_results$",
full.names =T)
deg_files
##读取一个差异基因分析结果
deg_results = read.table(file = deg_files[1]) %>%
rownames_to_column(var = "gene_id")
deg_results
##使用for循环结构批量读取差异基因分析结果
for(i in 2:length(deg_files)) {
d =read.table(file = deg_files[i])%>%
rownames_to_column(var = "gene_id")
deg_results = bind_rows(deg_results,d)
统计差异基因数目（如下图）
library(tidyverse)
select(deg_results,gene_id,sampleA,sampleB,logFC, logCPM, PValue ) %>%
mutate(direction = if_else(abs(logFC) < 1, 'NS',
if_else(logFC >= 1,'Up','Down'))) ->deg_results
##分组统计
filter(deg_results,direction != 'NS') %>%
group_by(sampleA,direction) %>%
summarise(count = n())
image.png
火山图代码
library(EnhancedVolcano)
EnhancedVolcano(
toptable = my_deg_gm,
lab = my_deg_gm$gene_id,
selectLab = key_gene,
x = "logFC",
y = "logCPM",
FCcutoff = 2,
pCutoff = 0.01
)
