如何使用（白嫖）他人文章数据进行批量转录组分析【上】

原创2023-03-02 21:10·积沐森林

零、下载数据

以这篇文献中的数据【The Plant Journel| 清华大学戚益军教授|水稻lncRNA多腺苷酸化的胁迫响应调控】为例，进行数据的搜索和下载。

  1. 在文献中搜索关键词【accession number】
    这个词一般就代表这篇文章数据上传到NCBI的数据编号，GSE是基因表达数据库，我们先找到这个文章的项目号（Bioproject）PRJNA340947，然后找到找SRA数据编号SRP083700，然后选择一个数据下载比如这个SRR4098209, 点Run→Data access，最后随便用什么下载上传到服务器就可以了。【推荐在服务器装一个aria2c 直接下载到服务器】


一、数据存储路径

/home/rice_blast/0_fastqc/ 建立一个文件夹，
记住你文件的存储路

二、效验数据完整性

  • 因为在你下载几十G的数据过程中可能会出现数据中断等情况，需要你对你的数据进行一下完整性验证，保证和当时上传的数据是一致的。这里就涉及一个概念MD5，可以用来验证数据完整性。
  • 因此你需要获得文章作者上传的原始MD5值，用于和你下载的数据比较，如果相同这个数据才可以放心使用。去这个网址【ENA Browser (ebi.ac.uk)】获得SRA数据的MD5值或者直接下载数据和其他相关信息也可以。
  • 还是以SRR4098209为例子，如Gif 图所示。在搜索框里搜索，然后点击ReadFile下的一个三角号，找到SRR_MD5或者Fastq_MD5，根据自己下载的数据选择，也可以直接搜索项目号PRJNA340947批量下载这个项目中所有SRR数据包括MD5在内的信息。
  • 最后点击Tsv格式导出即可。
MD5是一种常用的密码学哈希函数，它将任意长度的消息（Message）转换为128位的消息摘要（Message Digest）。
MD5散列值的应用非常广泛，其中之一就是校验文件完整性。
可以通过计算文件的MD5值，将其与预先生成的或者参考文档中给出的MD5值进行比较，以验证文件是否被篡改或损坏。
md5sum 是一种常用的校验工具，它可以通过对文件进行 MD5 哈希值计算，来判断文件内容是否完整无误。
在使用 md5sum 工具对文件进行校验时，会生成一个 32 位的哈希值，
如果计算出的哈希值与源文件的哈希值相同，则说明文件完整无误，否则说明文件已经被修改或损坏。



  1. 下面是具体的批量验证MD5的代码，需要一点点编程基础，就是将directory=后面的换成你储存在上面那个网站下载的SRR的文件绝对路径。
  2. sra_list.txt 为你下载的MD5信息，格式为每一行包含一个文件名和对应的md5值，用空格分隔。如下图

#在脚本中，您需要指定要验证的文件夹路径，
即将“/path/to/directory/”替换为您要验证的文件夹的路径。#以下是一个批量验证文件md5值的脚本，假设文件名和md5值已经记录在名为“sra_list.txt”的文件中，其中每一行包含一个文件名和对应的md5值，用空格分隔：
#!/bin/bash
# 指定要验证的文件夹路径
directory="/home/rice_blast/00_SRR/"
# 遍历文件夹下的所有文件并验证MD5值
for file in "$directory"/*; do
    if [ -f "$file" ]; then
        # 提取文件名
        filename=$(basename "$file")
        # 计算MD5值
        md5=$(md5sum "$file" | awk '{print $1}')
        # 读取MD5值文件中对应的MD5值
        expected_md5=$(grep "$filename" sra_list.txt | awk '{print $1}')
        # 比较MD5值
        if [ "$md5" == "$expected_md5" ]; then
            echo "$filename OK"
        else
            echo "$filename FAILED"
        fi
    fi
done
注意：文件夹为绝对路径
  • 以下为运行结果，表明我下载的数据和原始数据不匹配，可能无法使用。

一、压缩不用的SRR文件

for i in {05..07}
do
  pigz -9 -p 16 SRR91582$i
done
for i in {96..98}
do
  pigz -9 -p 16 SRR91581$i
done

二、解压SRR文件

  • SRR是NCBI Sequence Read Archive（SRA）数据库中使用的一种文件格式，用于存储高通量测序数据。SRR文件包含已压缩的测序数据和元数据信息，可以使用SRA工具包中的fastq-dump等工具进行解压和处理。
  • 这里使用fasterq-dump 因为其可以使用多个线程来解压数据，速度非常快（-e 16 使用16个线程）。




#conda安装软件
conda install -c bioconda sra-tools
# 指定SRR号
SRR_ARRAY=(SRR9158193 SRR9158194 SRR9158195 SRR9158199 SRR9158200 SRR9158201 SRR9158202 SRR9158203 SRR9158204 SRR9158208 SRR9158209 SRR9158210)
## 循环遍历 SRR 号数组，并在每个处理器上运行 fastq-dump 命令 
for SRR in ${SRR_ARRAY[@]}; 
do 
  fasterq-dump --gzip --split-3 -e 16 $SRR -O ~/rice_blast/00_fastq
done
#上面这个有点笨，可以使用下面的其他循环格式替换。

三、质控 fastqc

  • 以下为批量进行数据质量控制的脚本，既然你要做数据分析，必须得了解你的数据，因此将进行数据的质量控制，了解数据质量等各方面信息。通常使用的软件就是Fastqc。
#!/bin/bash

# 安装fastqc
conda install -c bioconda fastqc
# 设置fastq.gz文件所在的目录
fastq_dir="/home/rice_blast/00_fastq/"
# 设置FastQC结果输出目录
qc_dir="/home/rice_blast/00_fastq/"
# 遍历fastq.gz文件并运行FastQC
for file in ${fastq_dir}*.fastq.gz; 
do
  echo "Running FastQC on file: $file"
  fastqc -t 16 -o ${qc_dir} $file
done

四、去除接头和低质量

  • 质控没问题以后，数据中有些接头和低质量序列需要去除，fastp是一款处理序列性能和速度都非常好用的软件。


第一种方法。


#!/bin/bash
# 定义fastq.gz文件路径和输出目录
fastq_dir="/home/rice_blast/00_fastq/"
out_dir="/home/rice_blast/1_trim/"
# 定义一个文件列表，每个元素包含输入和输出文件名
file_list=(
  "SRR9158195_1.fastq.gz SRR9158195_2.fastq.gz SRR9158195.fastp_R1.fq.gz SRR9158195.fastp_R2.fq.gz SRR9158195.fastp.json SRR9158195.fastp.html"
  "SRR9158193_1.fastq.gz SRR9158193_2.fastq.gz SRR9158193.fastp_R1.fq SRR9158193.fastp_R2.fq SRR9158193.fastp.json SRR9158193.fastp.html"
  "SRR9158194_1.fastq.gz SRR9158194_2.fastq.gz SRR9158194.fastp_R1.fq SRR9158194.fastp_R2.fq SRR9158194.fastp.json SRR9158194.fastp.html"
  "SRR9158199_1.fastq.gz SRR9158199_2.fastq.gz SRR9158199.fastp_R1.fq SRR9158199.fastp_R2.fq SRR9158199.fastp.json SRR9158199.fastp.html"
  "SRR9158200_1.fastq.gz SRR9158200_2.fastq.gz SRR9158200.fastp_R1.fq SRR9158200.fastp_R2.fq SRR9158200.fastp.json SRR9158200.fastp.html"
  "SRR9158201_1.fastq.gz SRR9158201_2.fastq.gz SRR9158201.fastp_R1.fq SRR9158201.fastp_R2.fq SRR9158201.fastp.json SRR9158201.fastp.html"
  "SRR9158203_1.fastq.gz SRR9158203_2.fastq.gz SRR9158203.fastp_R1.fq SRR9158203.fastp_R2.fq SRR9158203.fastp.json SRR9158203.fastp.html"
  "SRR9158204_1.fastq.gz SRR9158204_2.fastq.gz SRR9158204.fastp_R1.fq SRR9158204.fastp_R2.fq SRR9158204.fastp.json SRR9158204.fastp.html"
  "SRR9158209_1.fastq.gz SRR9158209_2.fastq.gz SRR9158209.fastp_R1.fq SRR9158209.fastp_R2.fq SRR9158209.fastp.json SRR9158209.fastp.html"
  "SRR9158210_1.fastq.gz SRR9158210_2.fastq.gz SRR9158210.fastp_R1.fq SRR9158210.fastp_R2.fq SRR9158210.fastp.json SRR9158210.fastp.html")
# 检查输出目录是否存在，如果不存在则创建
if [ ! -d "$out_dir" ]; then
  mkdir "$out_dir"
fi
# 遍历文件列表并运行fastp
for f in "${file_list[@]}"; do
  arr=($f) # 将元素分割成输入和输出文件名
  # 运行fastp
  echo "Running fastp on files: ${arr[0]} and ${arr[1]}"
  fastp \
    --detect_adapter_for_pe \
    --in1 "$fastq_dir/${arr[0]}" \
    --in2 "$fastq_dir/${arr[1]}" \
    --out1 "$out_dir/${arr[2]}" \
    --out2 "$out_dir/${arr[3]}" \
    --json "$out_dir/${arr[4]}" \
    --html "$out_dir/${arr[5]}" \
    -w 4 \
    -q 30 \
    -u 20 \
    -c \
    -n 4 \
    -y \
    -f 1 \
    -z 7 \
    2> "$out_dir/${arr[0]}.fastp.log"
done
-c 是对overlap的区域进行纠错，所以只适用于pairend reads。-q 参数指定碱基质量的阈值 -u 参数指定一条序列中允许的低质量碱基的百分比，取值范围从0-100，如果序列中低质量碱基百分比超过了该阈值，这条序列就会被过滤掉；-n 参数指定一条序列中最多允许的N碱基的个数，如果超过这个数值，这条序列会被过滤掉。--length_required 指定最小长度，小于该长度的reads会被过滤掉; -length_limit 指定最大长度，大于该长度的reads也会被过滤掉，如果不希望进行长度过滤，可以添加 -w 参数指定并行的线程数。-j 指定json格式的报告文件， -h 指定html格式的报告文件。-z/ --compression 输出压缩格式。给定一个数字1-9，调整压缩比率和效率的平衡；-l 接一个长度值，小于这个长度reads被丢掉，默认是15，这个在处理非illumina测序数据时很有用。--detect_adapter_for_pe 检测双端adapter

第二种方法


#fastq.gz文件所在的目录
fastq_dir="/home/rice_blast/00_fastq/"
# 设置Fastp结果输出目录
out_dir="/home/rice_blast/1_trim/"
# 遍历fastq.gz文件并运行Fastp
for file in ${fastq_dir}*.fastq.gz;
do
  # 提取文件名（不含扩展名）
  filename=$(basename "$file" | cut -d. -f1)
  filename="${filename%.*}"
  # 运行Fastp
  echo "Running Fastp on file: $file"
  fastp --detect_adapter_for_pe \
--in1 "${fastq_dir}${filename}.fastq.gz" \
--in2 "${fastq_dir}${filename}.fastq.gz" \
--out1 "${out_dir}${filename}_1_fastp.fq.gz" \
--out2 "${out_dir}${filename}_2_fastp.fq.gz" \
--json "${out_dir}${filename}_fastp.json" \
--html "${out_dir}${filename}_fastp.html" \
-w 16 \
-q 30 \
-u 20 \
-c \
-y \
-f 1 \
-n 4 \
-z 7 \
2> "${out_dir}${filename}_fastp.log"
done

五、去除rRNA

  • RNA-Seq数据中含有大量的rRNA序列，如果不进行去除，会严重影响后续分析的准确性和效率。常见的去除rRNA的方法包括基于比对的方法和基于序列特征的方法。其中，基于比对的方法是指将RNA-Seq数据比对到rRNA数据库中，将比对到rRNA的序列去除掉。而基于序列特征的方法则是根据rRNA序列的特征，筛选出不包含rRNA序列的reads。

其中，基于比对的方法在rRNA数据库的选择、构建和比对参数的设置等方面需要谨慎，而基于序列特征的方法则相对简单，但在筛选条件的选择上也需要注意。在实际应用中，可以结合这两种方法进行rRNA去除，以提高去除效果和准确性。

  • 常用的去除rRNA的软件包括：bowtie、HISAT2、SortMeRNA、rRNASelector等。对于不同的样品和数据类型，选用合适的软件和参数进行rRNA去除是非常重要的。这里使用bowtie。

  • 这里使用的植物rRNA数据库是使用的华南农业大学陈程杰老师整理泛植物的rRNA序列库，提取的逻辑就是在【Archive (arb-silva.de)】rRNA数据中分别下载SSU和LSU rRNA数据，其中包括9个分类，提取其中的

  • Archaeplastida

  • 原始植物一类，就包含了所有植物的rRNA, python 脚本在下方，如果想提取其他物种只需要将
  • Eukaryota;Archaeplastida
  • 这俩个替换成自己所在的物种分类就可以了。
import os
# 下载文件
os.system("aria2c -j 20 https://www.arb-silva.de/fileadmin/silva_databases/release_132/Exports/SILVA_1328.1_SSUParc_tax_silva_trunc.fasta.gz")
os.system("aria2c -j 20 https://www.arb-silva.de/fileadmin/silva_databases/release_132/Exports/SILVA_138.1_LSUParc_tax_silva_trunc.fasta.gz")
# 解压缩文件
os.system("pigz -d *.gz")
# 读取fasta文件，提取Eukaryota;Archaeplastida的rRNA序列
with open("panPlant.rRNA.fa", "w") as out_file:
    with open("SILVA_132_SSUParc_tax_silva_trunc.fasta") as in_file:
        seq_id = ""
        sequence = ""
        for line in in_file:
            if line.startswith(">"):
                if "Eukaryota;Archaeplastida" in seq_id:
                    out_file.write(seq_id + "\n" + sequence + "\n")
                seq_id = line.strip()
                sequence = ""
            else:
                sequence += line.strip()
        if "Eukaryota;Archaeplastida" in seq_id:
            out_file.write(seq_id + "\n" + sequence + "\n")
            
    with open("SILVA_132_LSUParc_tax_silva_trunc.fasta") as in_file:
        seq_id = ""
        sequence = ""
        for line in in_file:
            if line.startswith(">"):
                if "Eukaryota;Archaeplastida" in seq_id:
                    out_file.write(seq_id + "\n" + sequence + "\n")
                seq_id = line.strip()
                sequence = ""
            else:
                sequence += line.strip()
        if "Eukaryota;Archaeplastida" in seq_id:
            out_file.write(seq_id + "\n" + sequence + "\n")
1. Amorphea：无定形生物（Amorphea）是一个由一些真核生物组成的分类群，其中包括动物、真菌和一些原生生物，如滑鞭毛虫和变形虫等。这些生物之间的联系在分子水平上已经得到了确认。


2. Archaeplastida：原始植物（Archaeplastida）是一类真核生物，包括植物、绿色藻类和红藻。这些生物共同具有叶绿素a和b、卡诺酸和细胞壁中的纤维素。


3. Cryptophyceae：隐藻目（Cryptophyceae）是一类真核生物，它们具有类似植物的叶绿体和藻类的核糖体，但其细胞壁由纤维素和多糖构成。


4. Discoba：圆盘虫亚门（Discoba）是一类原生生物，包括两类线虫：螺旋体和鞭毛虫，以及一个单细胞生物类群：吞噬者。


5. Excavata：Excavata是一类原生生物，其中许多成员具有线粒体，但其线粒体的形态和功能与其他真核生物不同。该类群包括许多寄生生物和真核生物中唯一能够进行磷酸化呼吸的线虫。


6. Haptophyta：钩藻门（Haptophyta）是一类真核生物，具有特殊的鞭毛和钩状突起，同时还具有可以形成有机碳酸盐的质体。


7. Incertae：不确定分类（Incertae sedis）指的是一些生物的分类地位未能明确确定，或者尚未被归入已知的分类群之中。


8. Picozoa：皮孔虫（Picozoa）是一类单细胞生物，其细胞大小通常在1-5微米之间，是海洋中的一类原生生物。


9. SAR：SAR超群是一个由三类真核生物组成的分类群，其中包括硅藻、裸藻、蓝藻、镰刀藻和变形虫等生物，其名称源于其中三类生物的英文名称的首字母缩写。
source /gss1/env/bowtie.env
#激活环境
source /gss1/env/bowtie.env
#fastp.gz文件所在的目录
rmRNA_dir="/home/rice_blast/1_trim/11_trim/"
# 设置去除rRNA结果输出目录
out_dir="/home/rice_blast/1_trim/"
# 遍历fastq.gz文件并运行Fastp
for file in ${rmRNA_dir}*.fq.gz;
do
  # 提取文件名（不含扩展名）
  filename=$(basename "$file" | cut -d. -f1)
  filename="${filename%.*}"
   bowtie -p 16 -v 1 \
   -x ~/2021/rRNA_database/panPlant_rRNA_bowtie_index \
   -1 "${rmRNA_dir}${filename}.fastp_R1.fq.gz" \
   -2 "${rmRNA_dir}${filename}.fastp_R2.fq.gz" \
   "${out_dir}${filename}.sam" \
   --un "${out_dir}${filename}_rmRNA.fq" ;
done

下期继续更新转录组分析【下】。

举报
评论 0
    热门：