siRNA数据分析步骤

Small RNA，指长度小于200bp的RNA序列，通常不能翻译为蛋白，但在基因表达的调控中起作用。
一般而言，小RNA可以分为miRNA, piRNA, siRNA, snRNA, snoRNA, rRNA, tRNA等。
sRNA-seq没有固定的pipeline，根据具体需求有不同的分析步骤。本文整理了siRNA的处理流程。

---

所需软件：fastqc fastp blast+ bowtie mirdeep2

以我的数据为例:
工作目录/pub4/zhangchao/Methylome_Soybean/sRNA-seq

原始数据存放在fastq/raw_data中：

5个样本，每个样本两个重复，单端测序

脚本文件位于scripts目录

1、原始数据质控

使用fastqc对原始数据质控，代码如下:
raw_fastqc.sh

#!/bin/sh
# rawdata fastqc
if [ -d ~/Methylome_Soybean/sRNA-seq/fastq/rwa_data/fastqc]
then
    echo "fastqc directory has existed."
else
    mkdir ~/Methylome_Soybean/sRNA-seq/fastq/rwa_data/fastqc
    echo "fastqc directory is newly created."
fi
cd ~/Methylome_Soybean/sRNA-seq/fastq/raw_data
ls *fastq | cut -d . -f 1 |while read id;
do
    fastqc -i ${id}.fastq -o ./fastqc
done

部分结果如下：

可见拿到的数据质量已经很高了，并且adapter也已经去除。
接下来主要对reads长度进行筛选

2、质量过滤

使用fastp进行筛选，代码如下：
fastp.sh

#!/bin/bash
# fastp
if [ -d ~/Methylome_Soybean/sRNA-seq/fastq/clean_data ]
then
    echo "clean_data directory has existed."
else
    mkdir ~/Methylome_Soybean/sRNA-seq/fastq/clean_data
    echo "clean_data directory is newly created."
fi
cd ~/Methylome_Soybean/sRNA-seq/fastq/raw_data
ls *fastq | cut -d . -f 1 |while read id;
do
    fastp -i ${id}.fastq \
    -o ../clean_data/${id}_clean.fastq \
    -q 20 -u 20 \   #reads中，超过20%的碱基质量低于20则切除
    -5 -3 \         #两端滑窗，默认4bp，质量低于20则切除
    -A \            #这个参数意思是跳过去接头的步骤，不加则自动去接头，-a 可以指定接头
    --length_required 18 \  #过滤掉小于18bp的reads
    --length_limit 30       #过滤掉大于30bp的reads
done

#对过滤后的clean reads进行二次质控
if [ -d ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/clean_fastqc]
then
    echo "clean_fastqc directory has existed."
else
    mkdir ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/clean_fastqc
    echo "clean_fastqc directory is newly created."
cd ~/Methylome_Soybean/sRNA-seq/fastq/clean_data
ls *clean_fastq | cut -d . -f 1 |while read id;
do
    fastqc -i ${id}.fastq -o ./clean_fastqc
done

3、过滤掉已知的其它小RNA

小RNA测序得到的reads中也包含了许多其它小RNA，如rRNA, tRNA, snRNA, miRNA等，需要将这些序列去除。
我们从Rfam和RNAcentral下载已注释的小RNA序列，用blast方法将它们筛除，剩余的序列进行下一步分析。

3.1 小RNA库下载

3.1.1 Rfam库

进入Rfam

其中共有4192个fasta文件

创建 ~/Methylation_soyeban/sRNA_seq/rfam目录，下载Rfam库的所有fasta文件并解压

mkdir ~/Methylation_soyeban/sRNA_seq/rfam
cd ~/Methylation_soyeban/sRNA_seq/rfam
curl -O ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/fasta_files/RF0[0001-4192].fa.gz
gunzip *gz

接下来，复制Rfam库的taxonomy信息，获得库中rRNA，tRNA，snRAN，miRNA等对应的RF文件编号

将rfam.txt文件放到目录~/Methylation_soyeban/sRNA_seq/library中，根据小RNA分成多个txt文件

cd ~/Methylation_soyeban/sRNA_seq/library
grep 'rRNA' rfam.txt >rfam.rrna.txt
grep 'tRNA' rfam.txt >rfam.trna.txt
grep 'sRNA' rfam.txt >rfam.srna.txt
grep 'snRNA' rfam.txt >rfam.snrna.txt
grep 'miRNA' rfam.txt >rfam.mirna.txt

根据分类的文件信息，从Rfam库中筛选出大豆各类小RNA的fasta文件

cd ~/Methylation_soyeban/sRNA_seq/rfam
rrna_no=`awk '{print $1}' ../library/rfam.rrna.txt|sed 's/$/&\.fa/g'`
cat $rrna_no|sed -n '/Glycine max/{p;n;p}' > ../library/rfam.gm.rrna.fa

trna_no=`awk '{print $1}' ../library/rfam.trna.txt|sed 's/$/&\.fa/g'`
cat $trna_no|sed -n '/Glycine max/{p;n;p}' > ../library/rfam.gm.trna.fa

srna_no=`awk '{print $1}' ../library/rfam.srna.txt|sed 's/$/&\.fa/g'`
cat $srna_no|sed -n '/Glycine max/{p;n;p}' > ../library/rfam.gm.srna.fa

snrna_no=`awk '{print $1}' ../library/rfam.snrna.txt|sed 's/$/&\.fa/g'`
cat $snrna_no|sed -n '/Glycine max/{p;n;p}' > ../library/rfam.gm.snrna.fa

mirna_no=`awk '{print $1}' ../library/rfam.mirna.txt|sed 's/$/&\.fa/g'`
cat $mirna_no|sed -n '/Glycine max/{p;n;p}' > ../library/rfam.gm.mirna.fa

3.1.2 RNAcentral库

进入RNAcentral
在搜索框输入物种，在左侧边栏依次点选不同类型的小RNA并下载fasta文件到library目录

本次需要rrna, trna, snrna snorna, mirna, premirna

3.1.3 合并fasta文件并建库

cd ~/Methylation_soyeban/sRNA_seq/library
cat rfam.rrna.fa central.rrna.fa > gm.rrna.fa
cat rfam.trna.fa central.trna.fa > gm.trna.fa
cat rfam.srna.fa rfam.snrna.fa central.snrna.fa central.snorna.fa > gm.snrna.fa
cat rfam.mirna.fa central.mirna.fa central.premirna.fa > gm.snrna.fa
cd ~/Methylation_soyeban/sRNA_seq/blastdb
makeblastdb -in ../library/gm.rrna.fa -dbtype nucl -title gm.rrna -out gm.rrna
makeblastdb -in ../library/gm.trna.fa -dbtype nucl -title gm.trna -out gm.trna
makeblastdb -in ../library/gm.snrna.fa -dbtype nucl -title gm.snrna -out gm.snrna
makeblastdb -in ../library/gm.mirna.fa -dbtype nucl -title gm.mirna -out gm.mirna

3.2 reads简并

测序结果中含有许多相同的小RNA序列，使用mirdeep2将相同的reads简并到一条，可以提高blast效率

mirdeep2的collapse_reads_md.pl可以根据输入的fasta文件和给定的3字符代号将重复reads计数，为了方便将文件名作修改为s11.fastq, s12.fastq, s21.fastq…

#!/bin/bash
cd ~/Methylome_Soybean/sRNA-seq/fastq/clean_data
ls *.fastq | cut -d . -f 1|while read id;
do
    echo ${id}
    awk '{if(NR%4 == 2){print">"$0"\n"$0}}' ${id}.fastq > ${id}.fasta
    collapse_reads_md.pl ${id}.fasta ${id} > ${id}_collapsed.fasta
done

得到结果如下图所示，表示这条序列有333106次重复

3.3 blastn比对，获得fasta文件

用blast-short进行序列比对，挑选出没有比对到小RNA库的序列作为putative siRNAs

#!/bin/bash
if [ -d ~/Methylome_Soybean/sRNA-seq/blastresult ]
then
echo "the directory for blast_result has existed."
else
mkdir ~/Methylome_Soybean/sRNA-seq/blastresult
echo "the directory for blast_result is newly created."
fi
cd ~/Methylome_Soybean/sRNA-seq/fastq/clean_data
ls *collapsed.fasta |cut -d _ -f 1|while read id;
do
    mkdir ~/Methylome_Soybean/sRNA-seq/blastresult/${id}_blast
    cd ~/Methylome_Soybean/sRNA-seq/blastresult/${id}_blast

    blastn -task blastn-short -query ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}_collapsed.fasta -out ${id}_snrna.blast -db ~/Methylome_Soybean/sRNA-seq/blastdb/gm.snrna -outfmt 6 -evalue 0.01 -num_threads 8
    blastn -task blastn-short -query ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}_collapsed.fasta -out ${id}_mirna.blast -db ~/Methylome_Soybean/sRNA-seq/blastdb/gm.mirna -outfmt 6 -evalue 0.01 -num_threads 8
    blastn -task blastn-short -query ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}_collapsed.fasta -out ${id}_rrna.blast -db ~/Methylome_Soybean/sRNA-seq/blastdb/gm.rrna -outfmt 6 -evalue 0.01 -num_threads 8
    blastn -task blastn-short -query ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}_collapsed.fasta -out ${id}_trna.blast -db ~/Methylome_Soybean/sRNA-seq/blastdb/gm.trna -outfmt 6 -evalue 0.01 -num_threads 8
    echo "blast on ${id} against gm rRNA/tRNA/snRNA/miRNA completed."
  #根据剩余序列ID提取fasta文件
    cat *.blast | awk '{print $1}' | sort | uniq > annotated.txt
    sed -n '/^>/p' ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}_collapsed.fasta | sed 's/^>//g' | sort > all_seq.txt
    sort all_seq.txt annotated.txt annotated.txt | uniq -u > unannotated.txt
    perl ~/Methylome_Soybean/sRNA-seq/scripts/getseq.perl unannotated.txt ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}._collapsed.fasta > ${id}.filtered.fasta
  #找回被简并的数据
    grep -v '>' ${id}.filtered.fasta > ${id}.filtered.list
    perl /pub4/zhangchao/Methylome_Soybean/sRNA-seq/scripts/getseq.pl ${id}.filtered.list ~/Methylome_Soybean/sRNA-seq/fastq/clean_data/${id}.fasta > ${id}.filtered.all.fasta
  #计算剩余小RNA的长度分布
    seqkit fx2tab -l -n ${id}.filtered.fasta |awk '{print $2}'|sort |uniq -c > ${id}.filtered.dis
    seqkit fx2tab -l -n ${id}.filtered.all.fasta |awk '{print $2}'|sort |uniq -c > ${id}.filtered.all.dis
done

上述getseq.perl脚本

#perl getseq.pl  Osa.list Osa.port > Osa.fa
#!usr/bin/perl -w
open ID, "$ARGV[0]" || die $!;
while (<ID>){
 chomp;
 $hash{$_}=1;
}
open SEQ, "$ARGV[1]" || die $!;
while(defined($line = <SEQ>)){
 if($line =~ />/){
   $key = (split /\s/,$line)[0];
   $key =~ s/>//g;
}
print "$line" if exists $hash{$key} == 1;
}

3.4 筛选出 24nt siRNA

#filter 24nt sRNA from fasta
#python orign.fasta filtered.fasta

import sys

fasta = open(sys.argv[1],"r")
filtered = open(sys.argv[2],'w')
i = 0
for line in fasta:
        if i%2 == 0:
                seqID = line.strip("\n")
        elif i%2 == 1:
                sequence = line.strip("\n")
                if len(sequence) == 24:
                        qual=line.strip("\n")
                        output = seqID + "\n" + sequence + "\n"
                        filtered.write(output)
        i += 1