JiaoYuan's Blog
之前那一篇文章主要讲的是一些知识与工具的用法，这次用六组数据进行分析，得到基因表达矩阵。
还是先从 NCBI 搜索数据，查看数据是双端测序的数据之后，多选条目，直接 Accession list
我使用的是下列数据
SRR25907783
SRR25907784
SRR25907785
SRR25907786
SRR25907787
SRR25907788
使用 sra-tools 批量下载数据，并且使用 nohup 把程序挂在后台下载
nohup prefetch --option-file acc_list.txt
再使用脚本进行拆分
#!/bin/bash
mkdir SRR
mv ./SRR*/*.sra ./SRR
cd SRR
nohup fastq-dump --gzip --split-3 SRR*.sra
这里千万不要在程序没有运行完成的时候就进行下一步操作，使用 top 可以查看后台进程
当程序运行完成后便可看到后台任务中已经没有前面运行的程序了
同一目录下的 nohup.out 文件中是后台进程的运行记录
拆分完成后就需要进行质量检测，使用通配符批量检测，并且将检测报告放在单独一个文件夹以便后面进行压缩
mkdir fastqc_report
fastqc SRR2590778*.fastq.gz -o ./fastqc_report
完成之后，将检测报告进行合并，以便查看
multiqc ./fastqc_report
依据检测报告对序列进行过滤，参数之前已经讲过，这里数据比较多，写一个 bash 脚本的 for 循环：
#!/bin/bash
for i in {3..8}
do
trimmomatic PE -threads 1 -phred33 SRR2590778${i}_1.fastq.gz SRR2590778${i}_2.fastq.gz -summary oryza_sativa_${i}.summary -baseout SRR2590778${i}.fastq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:15 MINLEN:36
done
这里是对 SRR25907783 到 SRR25907788 的数据进行处理，根据自己数据的情况将这里进行修改即可
下载水稻的参考基因组和注释文件进行 hisat2 比对
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.57.gff3.gz
把下载的文件解压后重命名，方便使用
gzip -d Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz
mv Oryza_sativa.IRGSP-1.0.dna.toplevel.fa oryza_sativa.fa
gzip -d Oryza_sativa.IRGSP-1.0.57.gff3.gz
mv Oryza_sativa.IRGSP-1.0.57.gff3 oryza_sativa.gff3
先构建 hisat2 索引
hisat2-build oryza_sativa.fa oryza_sativa
构建完成之后开始将 reads 比对到参考基因组，还是写成脚本进行
#!/bin/bash
for i in {3..8}
do
hisat2 -x hisat2_index/oryza_sativa -p 5 -1 SRR2590778${i}_1P.fastq.gz -2 SRR2590778${i}_2P.fastq.gz -S SRR2590778_${i}.sam
done
还是写成 bash 脚本一键进行，这里我将每次排序的结构按照 1-6 分别命名
#!/bin/bash
for i in {3..8}
do
samtools sort -n -@ 5 SRR2590778_${i}.sam -o SRR2590778_${i}
done
最后一步 featureCounts 生成基因计数表
#!/bin/bash
bam_files=(*.bam) # 将所有的 bam 文件作为变量，输入给 featureCounts
if [ ${#bam_files[@]} -gt 0 ]; then
featureCounts -T 5 -t exon -g Name -a Lolium_perenne.gff3 -o gene.counts -p "${bam_files[@]}"
else
echo "No BAM files found in the current directory."
fi
至此完成，拿到矩阵
个人习惯完成分析后把文件进行一下归类
最后将数据下载到本地进行下游分析
全流程分析代码可以在我的 GitHub 获取。