#该流程用于对SRR数据直接分析
#若使用fastq文件分析，则去除前面的对于SRR转换的过程
#此位置需要更改为自己的基因组信息
REF="/mnt/f/test/ref/B73_v4"
GTF="/mnt/f/test/ref/Zea_mays.B73_RefGen_v4.47.gtf"
threads=22
mkdir 00Rawdata
#下载原始数据 #axel为可以多线程下载，比wget下载速度快一些
axel -n 22 -o ./00Rawdata ftp.sra.ebi.ac.uk/vol1/fastq/SRR603/007/SRR6039707/SRR6039707_1.fastq.gz
axel -n 22 -o ./00Rawdata ftp.sra.ebi.ac.uk/vol1/fastq/SRR603/007/SRR6039707/SRR6039707_2.fastq.gz
#利用md5sum 验证数据完整性
#md5sum -c ./00Rawdata/file.md5new.txt > ./00Rawdata/md5check.txt
mkdir 01Cleandata 02Alignment 03Result
mkdir ballgown
#替换成自己的文件
for i in SRR6039725 SRR6039724
do
#转化SRR为fastq文件
#可以使用多线程
fasterq-dump --split-3 --threads ${threads} ${i} -O ./00Rawdata/
#使用单线程
#fastq-dump --split-3 ${i} -O ./00Rawdata/
#压缩文件
#单线程
#gzip ./00Rawdata/${i}_1.fastq
#gzip ./00Rawdata/${i}_2.fastq
#多线程,默认使用全部线程
pigz ./00Rawdata/${i}_1.fastq
pigz ./00Rawdata/${i}_2.fastq
# 去接头
# trimmomatic 可执行文件的位置
# 单端或双端测序，单端为SE， 双端为PE
# phred33 或者-phred64 指定碱基质量编码，如果未指定，则自动识别。Illumina 使用-phred33
# -threads 指定使用的线程数
# 指定处理后的数据所放置的位置，单端测序如前所示，双端测序如下： ./00Rawdata/${i}_1.fastq.gz ./00Rawdata/${i}_2.fastq.gz ./01Cleandata/${i}_1_clean.fastq.gz ./01Cleandata/${i}_1_unpaired.fastq.gz ./01Cleandata/${i}_2_clean.fastq.gz ./01Cleandata/${i}_2_unpaired.fastq.gz
#LEADING:20 \ 从起始端开始去除低质量的碱基，只要一个碱基低于阈值，就会切除该碱基，并调查下一个碱基
# TRAILING:20从末端移除低质量的碱基，只要一个碱基质量低于与之，就会切除该碱基，并调查下一个碱基
#SLIDINGWINDOW:4:15 当窗口内的平均质量低于阈值时，执行活动窗口提出，通过考虑多个碱基，单个质量差的碱基不会导致删除了在reads后期高质量的数据
#MINLEN:20 设置保留reads的最小长度
#java -jar /mnt/d/ubuntuSoftware/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 -threads $threads ./00Rawdata/${i}_1.fastq.gz ./00Rawdata/${i}_2.fastq.gz ./01Cleandata/${i}_1_clean.fastq.gz ./01Cleandata/${i}_1_unpaired.fastq.gz ./01Cleandata/${i}_2_clean.fastq.gz ./01Cleandata/${i}_2_unpaired.fastq.gz LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:20
#fastp使用方法
fastp --thread ${threads} -i ./00Rawdata/${i}_1.fastq.gz -I ./00Rawdata/${i}_2.fastq.gz -o ./01Cleandata/${i}_1_clean.fastq.gz -O ./01Cleandata/${i}_2_clean.fastq.gz -j ./01Cleandata/${i}.json -h ./01Cleandata/${i}.html
#比对到基因组
#--dta这个选项告诉 HISAT2 生成用于基因表达定量的数据（DTA，Documented Transcripts Alignment）,还有专门用于下游cufflink分析的--dta-cufflinks
#-x 选项后面是参考基因组的索引文件的路径。在这里，参考基因组的索引文件位于 ../ref/B73_v4 目录中。
#-p 选项指定了要使用的线程数，这里设置为 22，表示使用 22 个线程来执行比对任务，以提高处理速度。
#-1,-2 选项用于指定待比对的双端测序数据文件的路径。
# -S 选项指定了输出比对结果的 SAM 格式文件的路径和名称。SAM 文件是一种常见的比对结果文件格式，它包含了每个测序读取的比对信息。
hisat2 --dta -x ${REF} -p ${threads} -1 ./01Cleandata/${i}_1_clean.fastq.gz -2 ./01Cleandata/${i}_2_clean.fastq.gz -S ./02Alignment/${i}.sam
#将sam文件转化为bam文件
#view 是 samtools 的一个子命令，用于查看、转换和筛选 SAM/BAM 文件。
#-@ 指定线程数
#-S 选项指定了要处理的输入文件，即 SAM 格式的比对结果文件的路径。
# -b 选项告诉 samtools 将输入文件转换为 BAM 格式。BAM 是一种二进制格式，通常比 SAM 文件更紧凑，占用更少的存储空间，并且更容易处理。
#> 这部分代码将转换后的 BAM 文件的输出路径和文件名指定为 ./02Alignment/${i}.bam。
samtools view -@ ${threads} -S ./02Alignment/${i}.sam -b > ./02Alignment/${i}.bam
#为了节约空间，可以将不用的sam文件删除
sudo rm -r ./02Alignment/${i}.sam
#sort 是 samtools 的一个子命令，用于对 BAM 文件进行排序。
#-@指定线程数
# 这是要排序的输入文件，即未排序的 BAM 文件的路径。${i} 变量似乎用于迭代多个样本文件，因此 ${i} 会在循环中被实际的文件名替换。这里假设未排序的 BAM 文件是以 ${i}.bam 的形式存储在 02Alignment 目录中。
#-o 选项指定了排序后的 BAM 文件的输出路径和文件名。排序后的 BAM 文件将保存为 ./02Alignment/${i}_sorted.bam。${i} 变量也在输出文件名中使用，以保持结果文件的命名与输入文件相关。
samtools sort -@ ${threads} ./02Alignment/${i}.bam -o ./02Alignment/${i}_sorted.bam
#为了节约空间，可以将不用的bam文件删除
sudo rm -r ./02Alignment/${i}.bam
#这是 stringtie 工具的名称，它是一个用于转录组装和量化的常用工具。
#-p 选项指定了要使用的线程数，这里设置为 22，表示使用 22 个线程来执行装配和量化任务，以提高处理速度。
#-G 选项指定了参考基因组的注释文件的路径，这是一个 GTF 文件。stringtie 使用这个文件来辅助装配和注释转录本。
#-e 选项表示要生成输入 BAM 文件中每个样本的单个转录本集合。
#-B 选项表示生成一个包含每个样本的 Ballgown 格式输出。Ballgown 是一种用于分析和可视化 RNA-seq 数据的工具，可以用于分析基因表达和调查转录本的变化。
#-o 选项指定了装配和量化后的转录本 GTF 文件的输出路径和文件名。${i} 变量似乎用于迭代多个样本文件，因此 ${i} 会在循环中被实际的文件名替换。这里生成的 GTF 文件将保存在
# -A 选项用于生成基因丰度的 TSV 文件，记录每个样本中每个基因的丰度信息。这个文件的路径和文件名也包括 ${i} 变量，以保持结果文件的命名与输入文件相关。生成的 TSV 文件将保存在 ./ballgown/${i}/gene_abundances.tsv。
# 这是输入 BAM 文件的路径，即已经比对和排序的测序数据的文件。${i} 变量也在输入文件名中使用，以处理多个样本文件。
stringtie -p ${threads} -G ${GTF} -e -B -o ./ballgown/${i}/transcripts.gtf -A ./ballgown/${i}/gene_abundances.tsv ./02Alignment/${i}_sorted.bam
done
#计算reads数
featureCounts -T ${threads} -p -a ${GTF} -o ./FeatureCounts_gene.txt -t gene -g gene_id ./02Alignment/*_sorted.bam
# 提取stringtie产生的FPKM以及TPM值
Rscript -e "rm(list = ls()); setwd('./');
library(dplyr); library(stringr); library(data.table);
myresult = list();
myresult$FPKM = list.files('ballgown/') %>% lapply(function(x){
file = read.table(paste0('ballgown/', x, '/gene_abundances.tsv'), sep = '\t', header = T, quote = '', check.names = F) %>% select(c('Gene ID', 'FPKM'))
colnames(file)[2] = x; return(file)
}) %>% Reduce(x = ., function(x, y){ full_join(x, y, by = 'Gene ID') });
myresult$TPM = list.files('ballgown/') %>% lapply(function(x){
file = read.table(paste0('ballgown/', x, '/gene_abundances.tsv'), sep = '\t', header = T, quote = '', check.names = F) %>% select(c('Gene ID', 'TPM'))
colnames(file)[2] = x; return(file)
}) %>% Reduce(x = ., function(x, y){ full_join(x, y, by = 'Gene ID') });
openxlsx::write.xlsx(myresult, 'Counts.tpm_fpkm.xlsx', rowNames = F, overwrite = T);
"
原文地址：https://www.jianshu.com/p/e18dec9c3f60