转录组分析：预处理和质控

记录一下分析代码

# 生成日志存放目录
mkdir logFiles #日志需要添加时间，ts模块来源于moreutils：sudo apt install moreutils

1. NGS数据准备

# 有参转录组数据来源于PRJNA854395
Liu Y, Cai Y, Li Y, et al. Dynamic changes in the transcriptome landscape of Arabidopsis thaliana in response to cold stress[J]. Frontiers in Plant Science, 2022, 13: 983460.

# 创建原始数据存放目录，这里我已经用seqtk提前抽取原始数据5%的数据用于示例练习
mkdir -p 00_data && cd 00_data

du -ah --max-depth=1 # 查看一下文件大小
#############################################################
77M     ./A1_1.fq.gz
79M     ./A1_2.fq.gz
79M     ./A2_1.fq.gz
83M     ./A2_2.fq.gz
70M     ./A3_1.fq.gz
73M     ./A3_2.fq.gz
80M     ./B1_1.fq.gz
84M     ./B1_2.fq.gz
76M     ./B2_1.fq.gz
78M     ./B2_2.fq.gz
75M     ./B3_1.fq.gz
77M     ./B3_2.fq.gz
92M     ./C1_1.fq.gz
96M     ./C1_2.fq.gz
78M     ./C2_1.fq.gz
82M     ./C2_2.fq.gz
72M     ./C3_1.fq.gz
76M     ./C3_2.fq.gz
70M     ./D1_1.fq.gz
73M     ./D1_2.fq.gz
78M     ./D2_1.fq.gz
80M     ./D2_2.fq.gz
85M     ./D3_1.fq.gz
89M     ./D3_2.fq.gz
74M     ./E1_1.fq.gz
77M     ./E1_2.fq.gz
81M     ./E2_1.fq.gz
84M     ./E2_2.fq.gz
69M     ./E3_1.fq.gz
72M     ./E3_2.fq.gz
2.3G    .
#############################################################

# 手动处理生成样本分组信息sample_list.txt，放在根目录下
A   A1
A   A2
A   A3
B   B1
B   B2
B   B3
C   C1
C   C2
C   C3
D   D1
D   D2
D   D3
E   E1
E   E2
E   E3

2. 参考基因组准备

# 回到上级目录，创建参考基因组存放目录
mkdir -p 01_gem && cd 01_gem

# 下载和解压拟南芥基因组文件，来源ensembl
https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-59/fasta/arabidopsis_thaliana/dna/

wget -O TAIR10.fa.gz https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-59/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz

wget -O TAIR10.gff3.gz https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-59/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.59.gff3.gz

pigz -k -d *.gz -p 16

# 使用pyGFF.py将GFF转换为GTF和BED，同时提取出一些信息
python $HOME/calculate/workflow/rna_ref/pyGFF.py TAIR10.gff3 --src ensembl --info --bed

提取出的信息表格TAIR10_feature_info.csv

# 提取出基因信息
awk -F ',' 'NR==1 || $3 == "-"' TAIR10_feature_info.csv > TAIR10_genes.csv

提取出的信息表格TAIR10_genes.csv

对于模式植物，可以从多个方面补充一些注释信息，比如GO和KEGG

# 拟南芥GO注释，来源于GO官网
wget -O ath.gaf.gz https://current.geneontology.org/annotations/tair.gaf.gz

pigz -k -d ath.gaf.gz # 解压

awk -F '\t' '!/^!/ && $2 == $10 {print $2 "\t" $4 "\t" $5}' ath.gaf > ath.GO.txt 
# 提取非!行中第二列基因ID和第十列基因Name相同的行，并取出这些行的第二列、第四列、第五列：基因ID/从属关系/GO编号

ath.GO.txt

aria2c -s 16 -x 16 https://purl.obolibrary.org/obo/go/go-basic.obo 
#下载go-basic.obo文件

python $HOME/calculate/workflow/rna_ref/get_go_term.py 
# 借助goatools解析go-basic.obo文件，得到GO层级关系GO_Term_hierarchy.txt
# pip install goatools -i https://pypi.tuna.tsinghua.edu.cn/simple

GO_Term_hierarchy.txt

# 拟南芥KEGG注释，来源于KEGG官网
python $HOME/calculate/workflow/rna_ref/get_org_kegg_info.py --ath
# 得到ath.KEGG.Pathway.txt

ath.KEGG.Pathway.txt

到这里，原始数据及基因组基本准备完毕，可以进入后续质控

3. 数据质控及过滤

# 回到上级目录，创建质控目录
mkdir -p 02_qc

mkdir 02_qc/rawFq && mkdir 02_qc/cleanFq

# 创建软链接，如果无法创建软链接，直接复制文件过去
ln -s $PWD/00_data/*.fq.gz $PWD/02_qc/rawFq

# 进行原始数据的fastqc质控，fastqc版本：0.12.1
fastqc -t 8 02_qc/rawFq/*.fq.gz -o 02_qc/rawFq 2>&1 | ts '%Y-%m-%d %H:%M:%S' | tee $PWD/logFiles/fastqc_rawFq.log

# multiqc合并fastqc质控报告，multqic版本：1.25.1
multiqc 02_qc/rawFq -o $PWD/02_qc -n rawFq_report 2>&1 | ts '%Y-%m-%d %H:%M:%S' | tee $PWD/logFiles/multiqc_rawFq.log

multiqc报告

# seqkit统计基本指标，seqkit版本：2.8.2
seqkit stats -a -j 8 02_qc/rawFq/*.fq.gz > 02_qc/rawFq_basic_stats.txt

# 转换成表格
python $HOME/calculate/workflow/rna_ref/fqstats2excel.py 02_qc/rawFq_basic_stats.txt
#此脚本需要pandas(2.2.3)/openpyxl(3.1.5)/xlsxwriter模块

rawFq_basic_stats.xlsx

# fastp过滤，fastp版本：0.23.4
for file in $PWD/02_qc/rawFq/*_1.fq.gz; do
    # 提取样本名称
    sample=$(basename "$file" _1.fq.gz)
    # 运行fastp过滤双端文件
    fastp \
    -i "$PWD/02_qc/rawFq/${sample}_1.fq.gz" \
    -I "$PWD/02_qc/rawFq/${sample}_2.fq.gz" \
    -o "$PWD/02_qc/cleanFq/${sample}_1.clean.fq.gz" \
    -O "$PWD/02_qc/cleanFq/${sample}_2.clean.fq.gz" \
    -j "$PWD/02_qc/cleanFq/${sample}.filtered.fastp.json" \
    -h "$PWD/02_qc/cleanFq/${sample}.filtered.fastp.html" \
    -w 16 \
    -q 20 \
    -l 50 \
    -u 50
done 2>&1 | ts '%Y-%m-%d %H:%M:%S' | tee $PWD/logFiles/trim_fastp.log

# 统计过滤完数据信息，并转换成表格
python $HOME/calculate/workflow/rna_ref/fastp2stats.py $PWD/02_qc/cleanFq -o $PWD/02_qc -x cleanFq_basic_stats.xlsx

cleanFq_basic_stats.xlsx

使用到的脚本参考

https://gitee.com/weihao21/pygff/blob/master/pyGFF.py
https://gitee.com/weihao21/pipeline/tree/master/rna_ref