第一部分：脚本的初始设置与参数解析
#!/bin/bash
# 检查输入参数
if [ "$#" -lt 1 ]; then
echo "Usage: $0 -f <sample_file> -d <data_directory> -s <script_directory> -g <group_file>"
exit 1
fi
# 使用 R 语言的 argparser 包解析命令行参数
Rscript parse_args.R "$@" > args.sh
# 加载解析出的参数
if [ ! -f args.sh ]; then
echo "Error: args.sh not found, something went wrong in the argument parsing."
exit 1
fi
# 加载解析出的参数
source ~/yiyaoran/rna-seq/my_rna_seq/1516-ran/args.sh
# 打印所有变量以进行调试
echo "SAMPLE_FILE: $SAMPLE_FILE"
echo "DATADIR: $DATADIR"
echo "SCRIPTDIR: $SCRIPTDIR"
echo "GROUP_FILE: $GROUP_FILE"
# 设置脚本在出现错误时停止
set -e
讲解：
- 参数检查：脚本首先检查是否传递了足够的参数，至少需要1个参数。如果不满足条件，则打印使用方法并退出。
- 参数解析：使用
Rscript调用R脚本parse_args.R来解析命令行参数，并将结果保存为args.sh文件。 - 加载参数：如果解析成功，脚本会加载
args.sh中的变量，这些变量将用于后续分析（如SAMPLE_FILE,DATADIR等）。 - 调试输出：脚本通过
echo打印加载的参数值，便于调试和确认输入是否正确。 - 设置错误中断：
set -e命令确保在出现任何错误时，脚本立即停止执行。
第二部分：FastQC 质量控制
# 确保 FastQC 输出目录存在
#OUTPUT_DIR="./1.fastqc_results"
#mkdir -p $OUTPUT_DIR
# 读取样本文件并对每个样本进行 FastQC 分析
#while IFS= read -r sample; do
# if [[ -f "$sample" ]]; then
# echo "Processing $sample with FastQC..."
# fastqc "$sample" -o "$OUTPUT_DIR"
# else
# echo "File not found: $sample"
# fi
#done < "$SAMPLE_FILE"
#echo "FastQC analysis complete. Results are stored in $OUTPUT_DIR."
第三部分：数据质控处理（fastp）
# 设置工作目录和相关路径
workdir=$(pwd)
mkdir -p "$workdir/2.data_qc"
cd "$workdir/2.data_qc"
# 提取样本名称列表
declare -a samples=()
while IFS= read -r sample; do
# 从路径中提取样本名，假设样本名格式为 {sample}_R1.fq.gz 或 {sample}_R2.fq.gz
basename=$(basename "$sample")
sample_name="${basename%_*}"
if [[ ! " ${samples[*]} " =~ " ${sample_name} " ]]; then
samples+=("$sample_name")
fi
done < "$SAMPLE_FILE"
# 运行 fastp 进行质控处理
for i in "${samples[@]}"; do
echo "RUN CMD: fastp --thread 4 --qualified_quality_phred 10 \
--unqualified_percent_limit 50 \
--n_base_limit 10 \
-i $DATADIR/${i}_R1.fq.gz \
-I $DATADIR/${i}_R2.fq.gz \
-o ${i}_1.clean.fq.gz \
-O ${i}_2.clean.fq.gz \
--adapter_fasta $DATADIR/illumina_multiplex.fa \
-h ${i}.html -j ${i}.json"
fastp --thread 4 \ # 使用4个线程并行处理
--qualified_quality_phred 10 \ # 设置合格碱基的质量值阈值为10，低于该值的碱基视为不合格
--unqualified_percent_limit 50 \ # 如果一条read中超过50%的碱基不合格，则该read会被过滤掉
--n_base_limit 10 \ # 如果一条read中N碱基的数量超过10个，则该read会被过滤掉
-i $DATADIR/${i}_R1.fq.gz \ # 输入文件1（R1方向的测序数据，gzip压缩格式）
-I $DATADIR/${i}_R2.fq.gz \ # 输入文件2（R2方向的测序数据，gzip压缩格式）
-o ${i}_1.clean.fq.gz \ # 输出清理后的文件1（R1方向）
-O ${i}_2.clean.fq.gz \ # 输出清理后的文件2（R2方向）
--detect_adapter_for_pe \ # 自动检测并去除双端测序中的接头序列
-h ${i}.html \ # 生成HTML格式的质控报告
-j ${i}.json # 生成JSON格式的质控报告
done
# 质控数据统计汇总：
python $SCRIPTDIR/qc_stat.py -d "$workdir/2.data_qc/" -o "$workdir/2.data_qc/" -p all_sample_qc
#csvlook all_sample_qc.tsv
| SampleID | rawDataReads | cleanDataReads | rawDataBase | cleanDataBase | Q20 | Q30 | GC |
| ----------- | ------------ | -------------- | ------------- | ------------- | ----- | ----- | ----- |
| 1516-714T-2 | 46,811,722 | 46,810,452 | 7,021,758,300 | 6,948,101,453 | 98.55 | 95.95 | 55.32 |
| 1516-714W-1 | 46,263,762 | 46,262,678 | 6,939,564,300 | 6,875,725,469 | 98.36 | 95.39 | 55.27 |
| 1516-715T-4 | 41,997,012 | 41,995,142 | 6,299,551,800 | 6,235,252,339 | 98.58 | 96.02 | 55.21 |
| 1516-715W-3 | 41,437,680 | 41,436,064 | 6,215,652,000 | 6,149,566,109 | 98.33 | 95.37 | 56.29 |
echo "Data QC with fastp complete. Results are stored in $workdir/2.data_qc."
- 提取样本名称：
- 使用
while循环读取样本文件（$SAMPLE_FILE），通过提取文件名的前缀部分（假设格式为{sample}_R1.fq.gz）来生成样本名称列表。
- 使用
第四部分：Reads 与基因组进行比对（HISAT2 和 Samtools 处理）
# 定义参考基因组索引和基因组 BED 文件路径
REF_INDEX="/public/home/2022099/yiyaoran/rna-seq/my_rna_seq/ref/maize/Zea_mays.Zm-B73-REFERENCE-NAM-5.0.dna.toplevel"
GENE_BED="/public/home/2022099/yiyaoran/rna-seq/my_rna_seq/ref/maize/gene.bed"
cd $workdir # 回到工作目录
mkdir -p $workdir/3.map/hisat2 # 创建比对结果目录
cd $workdir/3.map/hisat2
# 运行 HISAT2 进行比对
for i in "${samples[@]}"; do
echo "RUN CMD: hisat2 -p 4 --rg-id=${i} --rg SM:${i} --rg LB:${i} --rg PL:ILLUMINA \
-x $REF_INDEX --dta \
-1 $workdir/2.data_qc/${i}_1.clean.fq.gz \
-2 $workdir/2.data_qc/${i}_2.clean.fq.gz \
-S ${i}.sam 2>${i}.summary"
##hisat2软件，链特异性文库设置
#--rna-strandness RF or FR（RF: read 1是反向链，read 2是正向链）
#非链特异性文库，不设置此参数
hisat2 -p 4 \ # 使用4个线程进行比对
--rg-id=${i} \ # 添加read组ID（read group ID），设置为当前样本ID
--rg SM:${i} \ # 添加样本名称（sample），设置为当前样本ID
--rg LB:${i} \ # 添加文库ID（library），设置为当前样本ID
--rg PL:ILLUMINA \ # 指定测序平台为Illumina
-x $REF_INDEX \ # 指定参考基因组的索引文件
--dta \ # 为转录组分析（downstream transcriptome analysis）优化比对
-1 $workdir/2.data_qc/${i}_1.clean.fq.gz \ # 输入清理后的R1方向的reads
-2 $workdir/2.data_qc/${i}_2.clean.fq.gz \ # 输入清理后的R2方向的reads
-S ${i}.sam \ # 输出SAM格式的比对结果文件
2>${i}.summary # 将比对的日志信息输出到summary文件中
done
#cat 1516-714T-2.summary
#23405226 reads; of these:
# 23405226 (100.00%) were paired; of these:
# 1221258 (5.22%) aligned concordantly 0 times
# 20874507 (89.19%) aligned concordantly exactly 1 time
# 1309461 (5.59%) aligned concordantly >1 times
# ----
# 1221258 pairs aligned concordantly 0 times; of these:
# 84254 (6.90%) aligned discordantly 1 time
# ----
# 1137004 pairs aligned 0 times concordantly or discordantly; of these:
# 2274008 mates make up the pairs; of these:
# 1292738 (56.85%) aligned 0 times
# 895745 (39.39%) aligned exactly 1 time
# 85525 (3.76%) aligned >1 times
#97.24% overall alignment rate
#- **总共有 23405226 条 reads**：
# - 其中，**100.00%**（即 23405226 条 reads）都是成对的 paired-end reads；
#- 在这些成对的 reads 中：
# - **5.22%**（1221258 对）完全没有比对上；
# - **89.19%**（20874507 对）正好比对上1次；
# - **5.59%**（1309461 对）比对上多于1次；
#- **1221258 对 reads** 没有成对的比对上（即 concordantly 0 次），其中：
# - **6.90%**（84254 对）在不一致的情况下比对上了1次（即 discordantly 1 次）；
#- **1137004 对 reads** 没有成对地比对上，也没有在不一致的情况下比对上，其中：
# - 这些未比对上的 reads 共有 **2274008 个单端 mates**（构成这些配对），其中：
# - **56.85%**（1292738 条 reads）没有比对上；
# - **39.39%**（895745 条 reads）比对上了1次；
# - **3.76%**（85525 条 reads）比对上了多于1次。
#- **总体比对率为 97.24%**：这是指所有 reads 中成功比对上的比例，说明数据的比对质量较高。
# sam 转换为 bam 并排序
for i in "${samples[@]}"; do
echo "RUN CMD: samtools sort --threads 4 -m 3G -o ${i}.bam ${i}.sam"
samtools sort --threads 4 -m 3G -o ${i}.bam ${i}.sam
done
# bam index 建立索引
for i in "${samples[@]}"; do
echo "RUN CMD: samtools index ${i}.bam"
samtools index ${i}.bam
done
第五部分：对 Map 结果进行 QC 分析
主要使用RSeQC的·geneBody_coverage.py工具进行
geneBody_coverage.py 是 RSeQC 工具包中的一个脚本，主要用于评估 RNA-seq 数据中基因覆盖的均匀性。它帮助研究者检测 RNA 样本是否在转录本的 5’ 端到 3’ 端具有均匀的覆盖，从而评估样本 RNA 是否发生了降解。这种分析非常重要，因为 RNA 降解会导致数据偏差，影响下游基因表达分析的准确性。
主要功能：
geneBody_coverage.py 通过对参考基因组的注释文件（通常是 GTF 格式）和 RNA-seq 比对结果（如 BAM 文件）的分析，生成一个覆盖曲线图，显示基因从 5’ 端到 3’ 端的覆盖情况。
常见参数：
-r：输入基因组注释文件（GTF 或 BED 格式），指定基因的区域。-i：输入 RNA-seq 的比对文件（BAM 格式），其中包含了比对到基因组上的 reads。-o：指定输出文件的前缀，脚本会生成相应的输出文件，包括覆盖情况的图表和数据。
输出结果：
- 覆盖曲线图（gene body coverage plot）：该图展示了基因组从 5’ 端到 3’ 端的覆盖分布。如果 RNA 没有显著降解，覆盖应该在基因全长上比较均匀。
- 数值文件：包含覆盖数据，方便进一步分析。
工作流程：
- 脚本首先将基因的转录本分成 100 个等长区间，从 5’ 端到 3’ 端计算每个区间的覆盖率。
- 然后，它汇总所有基因的覆盖情况，绘制覆盖曲线图。
- 如果 RNA 样本在处理过程中发生了降解，通常会看到覆盖从 5’ 端到 3’ 端逐渐下降的趋势。
使用示例：
geneBody_coverage.py -r Homo_sapiens.GRCh38.104.gtf -i sample1.bam -o sample1.genebody
cd $workdir # 回到工作目录
mkdir -p $workdir/3.map/map_QC # 创建QC结果目录
cd $workdir/3.map/map_QC
# 片段 inner size 分析，检测片段选择是否异常
for i in "${samples[@]}"; do
echo "RUN CMD: inner_distance.py -i $workdir/3.map/hisat2/${i}.bam -r $GENE_BED -o ${i}_inner_size"
inner_distance.py -i $workdir/3.map/hisat2/${i}.bam -r $GENE_BED -o ${i}_inner_size
done
# 基因覆盖度分析，检测RNA是否降解
for i in "${samples[@]}"; do
echo "RUN CMD: geneBody_coverage.py -r $GENE_BED -i $workdir/3.map/hisat2/${i}.bam -o ${i}.genebody"
geneBody_coverage.py -r $GENE_BED -i $workdir/3.map/hisat2/${i}.bam -o ${i}.genebody
done
echo "Read mapping and QC complete. Results are stored in $workdir/3.map."
基因覆盖情况图（Gene Body Coverage Plot）
解释：这张图显示了基因从5’端到3’端的覆盖率（Coverage）。图中横轴代表基因的百分比位置（从5’端到3’端），纵轴是覆盖率。
含义：这张图用来评估RNA样本是否发生了降解。理想情况下，基因的覆盖应该是均匀的，呈现出像钟形的分布。图中显示的覆盖从5’端到3’端逐渐升高，随后平稳，然后再到3’端逐渐下降。这种形状表明RNA样本的质量较好，没有显著降解。
mRNA片段插入大小分布图（mRNA Insert Size Distribution）
解释：这张图展示了插入片段大小（即文库中测序到的片段长度）的分布情况，横轴是插入大小（bp），纵轴是密度。
含义：片段大小分布用来评估文库构建的质量。正常的插入片段大小分布应符合预期，通常在一个范围内具有峰值。在这张图中，片段大小的均值为负值（-40.65 bp），这可能表明数据中有不对称的片段分布，可能是由于实验过程中某些步骤出错了，比如片段过短或异常的反转录反应。
第六部分：基因表达定量及结果展示
cd $workdir/ # 回到工作目录
mkdir -p 4.expression # 创建基因表达定量结果目录
cd 4.expression
# GTF 文件路径
GTF_FILE="/public/home/2022099/yiyaoran/rna-seq/my_rna_seq/ref/maize/Zea_mays.Zm-B73-REFERENCE-NAM-5.0.59.gtf"
GENE_LENGTH="/public/home/2022099/yiyaoran/rna-seq/my_rna_seq/ref/maize/gene_length.txt"
# 使用 htseq-count 进行基因表达定量
for i in "${samples[@]}"; do
echo "RUN CMD: htseq-count --format bam --order pos --mode union \
--stranded no --minaqual 10 --type exon \
--idattr gene_id $workdir/3.map/hisat2/${i}.bam $GTF_FILE > ${i}_gene.tsv"
htseq-count \
--format bam \ # 输入文件格式为 BAM
--order pos \ # 按照比对结果中的位置顺序进行计数
--mode union \ # 使用 "union" 模式进行计数，表示一个 read 覆盖多个 exon 时，按最大重叠区域进行计数
--stranded no \ # 指定为非链特异性测序数据
--minaqual 10 \ # 只计入测序质量分值大于或等于10的 reads
--type exon \ # 只对注释文件中的 exon 类型进行计数
--idattr gene_id \ # 使用基因ID（gene_id）作为基因的唯一标识
$workdir/3.map/hisat2/${i}.bam \ # 输入 BAM 文件（比对结果文件）
$GTF_FILE > ${i}_gene.tsv # GTF 文件作为注释，并将结果输出到对应的 gene TSV 文件中
#cat 1516-714T-2_gene.tsv |head -n 3
#Zm00001eb000010 149
#Zm00001eb000020 1729
#Zm00001eb000030 0
done
# 合并不同样品的基因表达定量结果
MERGE_CMD="python $SCRIPTDIR/merge_gene_count.py -p all_gene_count" # 初始化命令，指定合并的基因计数脚本和输出文件名
for i in "${samples[@]}"; do # 遍历样本数组中的每个样本
MERGE_CMD+=" -f ${i}_gene.tsv -l ${i}" # 为每个样本添加基因计数文件及对应标签参数
done
echo "RUN CMD: $MERGE_CMD"
eval $MERGE_CMD
# 基因表达定量结果展示
Rscript $SCRIPTDIR/fpkm_and_plot.R -i all_gene_count.tsv -l $GENE_LENGTH -o ./
echo "Gene expression quantification complete. Results are stored in $workdir/4.expression."
第七部分：基因差异表达分析（DESeq2）及可视化
cd $workdir/ # 回到工作目录
mkdir -p 5.deg # 创建差异表达基因（DEG）分析结果目录
cd 5.deg
# 使用 DESeq2 进行差异表达分析
Rscript ~/yiyaoran/rna-seq/my_rna_seq/scripts/deseq_analysis.r \
-i $workdir/4.expression/all_gene_count.tsv \
-g $GROUP_FILE \
-k $workdir/4.expression/all_gene_fpkm.tsv \
-r W --fdr 0.01 --fc 2 -p $(basename $GROUP_FILE .compare.txt)
# 绘制差异表达基因热图
Rscript ~/yiyaoran/rna-seq/my_rna_seq/scripts/heatmap.r \
-i $workdir/4.expression/all_gene_fpkm.tsv \
-l $(basename $GROUP_FILE .compare.txt).DEG.final.tsv \
-p $(basename $GROUP_FILE .compare.txt).deg_gene_heatmap \
-o ./
echo "Differential expression analysis complete. Results are stored in $workdir/5.deg."
使用R包DEseq2进行进行差异表达分析
library("DESeq2")
#==============================================
# 读取readcount矩阵文件 #
#==============================================
countdata <- read.csv("CountMatrix.csv",header = T,row.names = 1)
head(countdata, 10)
#生成分组信息
coldata <- read.csv("sample.csv",header = T)
#ID group
#1516-714W-1 W
#1516-715W-3 W
#1516-714T-2 T
#1516-715T-4 T
#关键步骤，生成Deseq对象
dds <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ group)
#==============================================
# 查看DESeqDataSet对象 #
#==============================================
dim(dds)
assay(dds)
assayNames(dds)
colSums(assay(dds))
rowRanges(dds)
colData(dds)
#过滤没有reads比对上的基因，所有reads数为零
nrow(dds)
dds <- dds[rowSums(counts(dds)) > 1 , ]
nrow(dds)
colSums(assay(dds))
barplot(colSums(assay(dds)),las=3,col=rep(rainbow(4),each=2))
#==============================================
# 差异表达计算 #
#==============================================
dep <- DESeq(dds)
res <- results(dep)
res
write.csv(x = res,file = "des.csv")
#筛选出差异表达基因，log2foldchange >=1或者<=-1,并且q值小于0.05的基因
res <- na.omit(res)
res0.05 <- res[(abs(res$log2FoldChange)>=1 & res$padj <=0.05), ]
nrow(res0.05)
#筛选出差异表达基因，log2foldchange >=2,并且q值小于0.01的基因
res0.01 <- res[(abs(res$log2FoldChange)>=1 & res$padj <=0.01), ]
nrow(res0.01)
write.csv(res0.05,file = "res0.05.csv")
#找出差异表达最显著的基因，按log2差异倍数排序
head(res0.05[order(res0.05$log2FoldChange,decreasing = TRUE), ])
head(res0.05[order(res0.05$log2FoldChange,decreasing = FALSE), ])
#==============================================
# 结果可视化 #
#==============================================
#MA图
plotMA(res, ylim = c(-10,10))
#火山图
m <- res
head(m)
m <- na.omit(m)
plot(m$log2FoldChange,m$padj)
plot(m$log2FoldChange,-1*log10(m$padj))
plot(m$log2FoldChange,-1*log10(m$padj),xlim = c(-10,10),ylim=c(0,100))
m <- transform(m,padj=-1*log10(m$padj))
down <- m[m$log2FoldChange<=-1,]
up <- m[m$log2FoldChange>=1,]
no <- m[m$log2FoldChange>-1 & m$log2FoldChange <1,]
plot(no$log2FoldChange,no$padj,xlim = c(-10,10),ylim=c(0,100),col="blue",pch=16,
cex=0.8,main = "Gene Expression",
xlab = "log2FoldChange",ylab="-log10(Qvalue)")
points(up$log2FoldChange,up$padj,col="red",pch=16,cex=0.8)
points(down$log2FoldChange,down$padj,col="green",pch=16,cex=0.8)
abline(v = c(-1,1),h = 2)
脚本讲解（未完待续）
deseq_analysis.r
# 加载 getopt 包，用于解析命令行参数
library(getopt)
第二部分：抽象样本分组信息
# 抽象样本分组矩阵函数
abstract_factors <- function(file) {
# 读取分组信息文件
f_mat <- read.table(file, header = TRUE, check.names = FALSE)
# 返回分组矩阵
return(f_mat)
}
第三部分：双条件差异表达分析函数
# 定义双条件差异表达分析函数
std_two_condition_comp_deseq2 <- function(f_mat, args, data) {
# 加载DESeq2库
library("DESeq2")
# 提取分组信息
condition <- f_mat[, 2]
print(condition)
# 确认分组是否为两个水平
level <- length(levels(as.factor(condition)))
if (level != 2) {
stop("分组的水平数量必须等于2")
}
# 创建差异表达分析对象
con <- as.factor(condition)
print(con)
# 创建DESeq2的数据集对象
colData <- data.frame(condition = f_mat[, 2])
rownames(colData) <- f_mat[, 1]
dds <- DESeqDataSetFromMatrix(countData = data,
colData = colData,
design = ~ condition)
# 重新设定参考组
dds$condition <- relevel(dds$condition, ref = opt$ref)
# 运行差异表达分析
dds <- DESeq(dds)
res <- results(dds, cooksCutoff = FALSE)
# 输出差异表达分析结果
print("差异表达分析结果：")
print(head(res))
# 过滤显著差异表达的基因
# 判断 fold change (FC) 的阈值：
# 如果传入的 opt$fc 参数小于等于 1，则将 fc_cutoff 设为 0；否则，将其设为 log2 转换后的 opt$fc 值
fc_cutoff <- ifelse(as.numeric(opt$fc) <= 1, 0, log2(as.numeric(opt$fc)))
# 判断差异表达基因（DGE）的条件：
# 1. res$padj 小于给定的 FDR 阈值（false discovery rate）
# 2. res$log2FoldChange 的绝对值大于 fc_cutoff（即变化倍数超过设定的阈值）
isDGE <- (res$padj < as.numeric(opt$fdr)) & (abs(res$log2FoldChange) > fc_cutoff)
# 标记显著性
# 将 res$padj 转换为向量，并存储在 significant 变量中
significant <- as.vector(res$padj)
# 将所有 NA 值替换为 "Normal"
significant[is.na(significant)] <- "Normal"
# 将符合显著性阈值（padj < fdr）并且 log2FoldChange 大于 fc_cutoff 的基因标记为 "Up"
significant[(res$padj < as.numeric(opt$fdr)) & (res$log2FoldChange > fc_cutoff)] <- "Up"
# 将符合显著性阈值（padj < fdr）并且 log2FoldChange 小于 -fc_cutoff 的基因标记为 "Down"
significant[(res$padj < as.numeric(opt$fdr)) & (res$log2FoldChange < -fc_cutoff)] <- "Down"
# 对于不显著（padj >= fdr）或 log2FoldChange 在 [-fc_cutoff, fc_cutoff] 之间的基因，标记为 "Normal"
significant[(res$padj >= as.numeric(opt$fdr)) | ((res$log2FoldChange <= fc_cutoff) & (res$log2FoldChange >= -fc_cutoff))] <- "Normal"
# 将 significant 转换为因子类型
significant <- as.factor(significant)
# 输出过滤后的差异表达基因数
print(paste("显著差异表达基因数量：", sum(isDGE), sep = ""))
# 返回差异表达基因数量
return(length(isDGE))
}
