1.建立相应目录
对新数据建立对应实验人员(xieruyu)、测序类型(RNA_seq)和日期(2022_07_28)的目录。
# 新建对应的目录
mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt scripts_log
2.检查数据完整性
(base) [customer@node01 2022_07_28]$ cat MD5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
3.在Linux服务器中对RNA_seq数据进行处理
vim新建RNA_seq_script将数据质控、清洗、比对、格式转换、排序、拼接和定量综合在一起。
#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法,当程序被执行时,能够载入bash的相关环境配置文件。
# Program
# This program is used for RNA-seq data analysis.
# History
# 2022/09/12 zexing First release
# 设置变量${dir}为常用目录
dir=/home/customer/lizexing/projects/xieruyu/RNA_seq/2022_07_28
# 对数据进行质控
fastqc -t 4 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz
# 利用for循环进行后续操作
for i in A1_FRAS220122137 A2_FRAS220122138 A3_FRAS220122139 B1_FRAS220122140 B2_FRAS220122141 B3_FRAS220122142 H1_FRAS220122143 H2_FRAS220122144 H3_FRAS220122145 I1_FRAS220122146 I2_FRAS220122147 I3_FRAS220122148 J1_FRAS220122149 J2_FRAS220122150 J3_FRAS220122151 K1_FRAS220122152 K2_FRAS220122153 K3_FRAS220122154 L1_FRAS220122155 L2_FRAS220122156 L3_FRAS220122157
do
# 对数据进行接头处理
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 -j 4 --paired \
${dir}/clean_data/${i}-1r_1.clean.fq.gz \
${dir}/clean_data/${i}-1r_2.clean.fq.gz
# 对数据进行比对
hisat2 -t -p 4 -x /home/customer/lizexing/references/NCBI/Gallus/hisat2_index/hisat2_index \
-1 ${dir}/clean_data/${i}-1r_1.clean.fq.gz \
-2 ${dir}/clean_data/${i}-1r_2.clean.fq.gz \
-S ${dir}/sam/${i}.sam
# 对数据进行格式转换
samtools view -@ 4 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam
# 对数据进行排序
samtools sort -@ 4 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam
# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 4 -G /home/customer/lizexing/references/NCBI/Gallus/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done
后台运行RNA_seq_script:
nohup bash RNA_seq_script > RNA_seq_script_log &
4.使用prepDE.py脚本提取read_counts数值
- 进入ballgown文件夹,将prepDE.py脚本拷贝至当前文件夹
cp /home/customer/lizexing/software/prepDE.py3 ./
- 退出当前conda环境
conda deactivate
- 使用python命令直接运行脚本
python prepDE.py3
运行结果中"gene_count_matrix.csv"即是DESeq2的输入文件。