一边学习，一边总结，一边分享！
由于微信改版，一直有同学反映。存在长时间接收不到公众号的推文。那么请跟随以下步骤，将小杜的生信筆記设置为星标，不错过每一条推文教程。
欢迎关注《小杜的生信笔记》！！
如何加入社群
小杜的生信笔记，仅有微信社群。
1. 微信群：付费社群。添加小杜好友，加友请知：加友须知！！，加入社群请查看小杜生笔记付费加友入群声明。
2. 小杜个人微信：若你有好的教程或想法，可添加小杜个人微信。值得注意的是，小杜个人微信并不支持免费咨询长时间咨询，但支持小问题2-3个免费咨询。
小杜微信：
知识星球：
2022年教程总汇
2023年教程总汇
1. 引言
大家好，我是小杜。
我们有1周未做更新了。最近，自己事情比较多，以及最近两周也在“水”一篇文章，今天刚刚投出去。因此，没时间来更新相关的教程。
若是很早就关注，我们的同学可能知道，我这边的更新都是看心情更新的，以及结合学习的情况来做更新，因为，自己一个人实在没有太多的精力和时间来的学习自己不涉及的内容。
因此，若是你“有”以及“愿意”分享相关的教程，可以投稿到我们这里。小杜也会给与适当的稿酬，一般在50-150之间。
最近，这一个月以来，基本是早上下地，上午上班搞东西。
今天的教程来自我们社群同学“immortals”的投稿。关于，samll RNA上游分析的整套流程，也主要是基于在自己的学习笔记进行记录，仅供大家参考。
致谢：在此感谢“immortals”同学的整理和投稿。此外，具投稿人说明也感谢
Bioinfor生信云公众号；以及，也感谢我们知识分享博主的分享。
本教程，我们会分几个小的章节免费提供给大家，大家可以进行订阅。
若是教程有错误的地方，请大家指出。
2. Small RNA教程
2.1 教程中涉及的流程
原始数据的质量控制使用 fastqc，使用trim_galore进行去除接头与低质量reads，并保留18-30nt reads。将数据 mirdeep2软件包中mapper.pl程序将clean reads从fq转成无冗余的fa格式。对各个长度sRNA数据的长度进行统计。
2.2 Small RNA分析中涉及到的相关数据库
将去除冗余后的数据，比对到 基因组上的数据进行miRNA分类注释。比对到·Rfam·、本物种miRNA的前体、本物种的repeat(重复序列)、本物种的mRNA、本物种的intron。（注释顺序·known mirna > ncrna > repeat > RNA > intron）。
下载Rfam数据库，去除其中的 miRNA，即为Rfam_rmMIR.fa(rRNA、tRNA、snRNA 和 snoRNA)。（也可生成一个Rfam_MIR.fa） 下载物种的成熟miRNA、前体miRNA。本次研究的是植物，可以在https://www.pmiren.com/下载植物数据或者www.plantsRNAs.org。 使用RepeatMasker进行重复序列注释，得到gff，再用bedtools进行提取。 根据 TBTools中插件getIntronGXF提取intron位置信息，参考https://mp.weixin.qq.com/s/KQBe4IlMpD2qQqJqwIp8aQ。 进行mRNA序列提取，根据genome.gff文件。 根据注释结果进行各首位碱基和各位碱基进行统计。
2.3. 定量
根据注释结果只选择 unkonw 与 intron的进行 novel miRNA进行预测及定量。
本文研究的是植物，采用 mirDP-P2进行novel miRNA预测。分别合并已知和新miRNA前体和成熟体序列，使用 mirdeep2中的quantifile.pl进行定量，RdDM分析时可以采用 ShortStack， 然后使用ShortStack软件将这些reads 映射到参考基因组 ，并只分析24-nt siRNA（mismatches=0,dicermin=24, dicermax=24）。
Shortstack 分析时，对不同前miRNA定量结果会在在一起，mirdeep2会将不同前体的miRNA结果分别定量。
2.4 靶基因预测
使用psRobot和targetfinder两款软件进行植物靶基因预测。
3. 软件的安装
安装所需的软件，在这里我们也可以创建一个新的环境，用于运行samll RNA。
mamba install -y trim-galore fastqc mirdeep2 bowtie samtools deeptools
4. samll RNA上游分析流程
4.1 过滤接头
# 对raw data进行 fastqc
ls raw_data/*.gz | while read line;
do fastqc -f fastq -o ./ $line;
done
# 保留 18nt~30nt 的序列，去除接头
trim_galore --length 18 --max_length 30 --stringency 3 --phred33 --cores 4 --dont_gzip -o clean_data2/ raw_data/Bud_1.fq.gz
# 过滤接头低质量reads后进行质控
ls clean_data2/*fq |while read line; do fastqc -f fastq -o clean_data2/ $line; done
4.2 合并冗余序列
# 合并冗余序列 从 id 部分反映每种序列出现的次数
# mirdeep2 软件包中mapper.pl程序
mapper.pl clean_data2/Bud_1_trimmed.fq -e -g B01 -h -m -s uniq_data/B01.fa &
ls clean_data2/*.fq | while read fq_file; do
# 提取文件名中的部分信息作为参数
base_name=$(basename "$fq_file" .fq)
sample_id=$(echo "$base_name" | cut -d'_' -f2)
# 构建 uniq_data 目录中的对应文件路径
uniq_fa="uniq_data/${sample_id}.fa"
#运行 mapper.pl 命令
mapper.pl "$fq_file" -e -g "$sample_id" -h -m -s "$uniq_fa" &
done
# 等待所有后台任务完成
4.3 长度统计
perl ../script/stat_srna.pl uniq_data/B01.fa B01
## 生成结果
(base) [yanghj@cln01 S1.reads_filter]$ head B01.readstat
sample total_read total_base uniq_read uniq_base
B01 23000873 528614567 8971764 207829072
(base) [yanghj@cln01 S1.reads_filter]$ head B01.len.total.txt
B01 18 481749
B01 19 575862
B01 20 876193
B01 21 2464188
B01 22 4121709
B01 23 2810818
B01 24 10073982
B01 25 524727
B01 26 381374
B01 27 186500
(base) [yanghj@cln01 S1.reads_filter]$ head B01.len.uniq.txt
B01 18 177427
B01 19 245442
B01 20 364507
B01 21 614174
B01 22 824131
B01 23 1564744
B01 24 4611942
B01 25 282459
B01 26 80698
B01 27 61581
脚本：stat_srna.pl代码。
#!/usr/bin/env perl
use strict;
use Bio::SeqIO;
die "perl $0 <srna.fasta> <out.prefix> \n" unless @ARGV ==2;
my $infile = shift;
my $outpre = shift;
my $total_read = 0;
my $total_base = 0;
my $uniq_read = 0;
my $uniq_base = 0;
my %len;
my $fa = Bio::SeqIO->new(-file => $infile, -format => 'fasta');
while(my $obj = $fa->next_seq()){
my $id = $obj->id;
my $seq = $obj->seq;
my $len = $obj->length;
## P11_879699_x139853
my @id = split /_/, $id ;
$id[2] =~ s/^x// ;
$total_read += $id[2];
$total_base += $id[2] * $len;
$uniq_read ++ ;
$uniq_base += $len;
$len{$len}{"uniq"} ++;
$len{$len}{"total"} += $id[2] ;
}
## output reads stats
open O1, ">$outpre.readstat" or die $!;
print O1 "sample\ttotal_read\ttotal_base\tuniq_read\tuniq_base\n";
print O1 "$outpre\t$total_read\t$total_base\t$uniq_read\t$uniq_base\n";
close O1;
## output length distribution
open OT, ">$outpre.len.total.txt\n";
open OU, ">$outpre.len.uniq.txt\n";
my @len = sort {$a<=>$b} keys %len;
my $min = $len[0];
my $max = $len[-1];
for(my $l = $min; $l <= $max ; $l++){
my $t = 0 ;
my $u = 0 ;
if(exists $len{$l}){
$t = $len{$l}{"total"} ;
$u = $len{$l}{"uniq"} ;
}
print OT "$outpre\t$l\t$t\n";
print OU "$outpre\t$l\t$u\n";
}
close OT;
close OU;
4.4 合并单个txt结果文件
将单个txt合并，绘制一张pdf
# 根据 *.len.total.txt 在R中自己绘图
cat *.len.uniq.txt > all.len.uniq.txt
Rscript script/draw_length_srna.R all.len.uniq.txt all.len.uniq.txt
cat *.len.total.txt > all.len.total.txt
Rscript script/draw_length_srna.R all.len.total.txt all.len.total.txt
脚本：draw_length_srna.R代码
#!/usr/bin/env Rscript
library(argparser, quietly=TRUE)# Create a parser
p <- arg_parser("Draw Length distribution for sRNAs")
# Add command line arguments
p <- add_argument(p, "length", help="input: length stat file, 3 column", type="character" )
p <- add_argument(p, "output", help="output prefix", type="character")
# Parse the command line arguments
argv <- parse_args(p)
input <- argv$length
output <- argv$output
#input <-"P11.len.total.txt"
#output <- "P11.len.total.txt"
len <- read.delim(input, header = F)
library(ggplot2)
pdf(paste(output, "length.pdf", sep = "."), width = 10, height = 7)
ggplot( len , aes(x = V2, y = V3)) +
geom_bar(
stat = 'identity',
aes(fill = V1),
position = 'dodge',
width = .5) +
theme_classic() +
ylab('count') +
xlab(' ') +
labs(fill = ' ') +
scale_x_continuous(breaks = c(min(len$V2) : max(len$V2)))
dev.off()
若我们的教程对你有所帮助，请点赞+收藏+转发，大家的支持是我们更新的动力！！
往期部分文章
1. 最全WGCNA教程（替换数据即可出全部结果与图形）
推荐大家购买最新的教程，若是已经购买以前WGNCA教程的同学，可以在对应教程留言，即可获得最新的教程。（注：此教程也仅基于自己理解，不仅局限于此，难免有不恰当地方，请结合自己需求，进行改动。）
2. 精美图形绘制教程
3. 转录组分析教程
4. 转录组下游分析
小杜的生信筆記 ，主要发表或收录生物信息学教程，以及基于R分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料!!
