一边学习，一边总结，一边分享！
由于微信改版，一直有同学反映。存在长时间接收不到公众号的推文。那么请跟随以下步骤，将小杜的生信筆記设置为星标，不错过每一条推文教程。
欢迎关注《小杜的生信笔记》！！
如何加入社群
小杜的生信笔记，仅有微信社群。
1. 微信群：付费社群。添加小杜好友，加友请知：加友须知！！，加入社群请查看小杜生笔记付费加友入群声明。
2. 小杜个人微信：若你有好的教程或想法，可添加小杜个人微信。值得注意的是，小杜个人微信并不支持免费咨询长时间咨询，但支持小问题2-3个免费咨询。
小杜微信：
知识星球：
2022年教程总汇
2023年教程总汇
5. miRNA进行分类注释
5.1 比对到参考基因组的reads保留，在进行以下的分类注释
# 安装比对软件(注意：若以安装，略过此步骤）
mamba install bowtie samtools deeptools
## reads比对基因组
bowtie -f \ # reads文件格式为fasta
-v 0 \ # 允许的错配数
-k 1 \ # 每条reads输出比对次数
--al P11.mapped.fa \ # 比对上的reads单独生成文件，进行后续分类注释
../S0.ref_prepare/03.genome/genome \ # bowtie index
../S1.reads_filter/uniq_data/P11.fa \ # 输入reads
-S P11.mapped.sam \ # 方便后期导入IGV查看，先要排序 及转化成bam
> P11.genome.bwt \ # 比对结果文件
2>P11.genome.log # 日志
# 将sam转化为bam，方便在IGV查看比对基因组情况文件
ls *mapped.sam | while read line; do
s=$(basename "$line" .fa)
samtools sort -O bam -@ 20 -o ${s}_sorted.bam ${s}.fa.sam &&\ # 进行排序
samtools index ${s}_sorted.bam && \ # 进行创建index
bamCoverage -b ${s}_sorted.bam -o ${s}.bw -p 15 --normalizeUsing RPKM --binSize 10 --extendReads 200 & # 转化成bw
done
wait
5.2 比对本物种miRNA前体数据库
5.2.1 比对本物种miRNA前体数据库
bowtie -f -v 0 \ # 允许的错配数
-p 5 \ # 比对线程数，根据自己情况设定
-a \ # 输出所有比对结果
--best --strata \ # 输出比对最好的结果
../S0.ref_prepare/01.known_miRNA/hairpin.ccr \ # database文件
./P11.mapped.fa \ # reads文件
> P11.miRNA.bwt 2>P11.miRNA.log
5.2.2 本物种成熟miRNA、前体的下载
https://www.pmiren.com
www.plantsRNAs.org
5.3 比对Rfam_rmMIR.fa
5.3.1 Rfam_rmMIR.fasta数据的制作
1. 下载Rfam数据
# 下载Rfam数据库中
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.4/database_files/family.txt.gz
wget -np -nH -r ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.4/fasta_files/
2. 将family文件中的type和家族信息进行提取
cut -f 1,2,19 $family > family.txt1
3. 整理Rfam中非miRNA家族数据,删除其中的miRNA家族
awk '$0 !~ /Gene; miRNA/' family.txt1 > family.txt.rmMIR
4. 获得Rfam_rmMIR.fasta
# 对家族列表中每一个家族fasta序列文件
# 在其序列id前添加家族编号，方便后续分析
cut -f 1 family.txt.rmMIR | while read line ;do \
zcat $fastadir/$line.fa.gz| sed "s/^>/>$line:/" \
done > Rfam_rmMIR.fasta
5.3.2 使用bowtie建立Rfam_rmMIR.fasta索引
bowtie-build Rfam_rmMIR.fasta Rfam_rmMIR.fasta
5.3.3 基于bowtie比对ncRNA数据库
bowtie -f -v 2 -p 1 -a --best --strata \ # 输出比对最好的结果
../S0.ref_prepare/02.Rfam/Rfam_rmMIR.fasta \ # bowtie index
./P11.mapped.fa > P11.ncRNA.bwt 2>P11.ncRNA.log
5.3.4 构建Rfam blast库
### 循环处理miRNA家族序列，添加家族编号并合并
cut -f 1 family.txt.MIR | while read line ;do \
zcat $fastadir/$line.fa.gz| sed "s/^>/>$line:/" \
done > Rfam_MIR.fasta
### 构建Rfam blast 库
makeblastdb -in Rfam_MIR.fasta -dbtype nucl
5.4 比对到本物种的重复序列
1. 获得repeat.fa文件
在这这里，作者提供了几种方式，获得repeat.fa文件，结合自己可以运行其中一种方法即可。
方法一：使用RepeatMasker快速比对生成repeat.gff文件
# RepeatMasker 软件自己安装难度有点高
# 没有repeat注释结果，可以自己注释之后生成repeat位置信息文件
RepeatMasker -pa 2 \ # 设置线程数2
-qq \ # 快速模式
-species all \ # 设置物种
-gff \ # 生成repeat位置的gff文件
-small \ # 生成的序列文件用小写字母代替repeat区序列
genome.fa >repeatmasker.log
# 将gff 1-based 转化为 bed 0-based
awk '{print $1"\t"$4-1"\t"$5}' repeat.gff > repeat.bed
使用bedtools进行提取。
bedtools getfasta -fi genome.fasta -bed repeat.bed -fo repeat.fa -name
方法二：基于singualarity，使用GenekTEtools进行重复序列注释
singualarity安装，需要管理员权限，如何没有请联系服务器管理员。或者使用
EDTA进行注释
# 构建数据库
singularity exec --env SINGULARITYENV_LANG=C --env SINGULARITYENV_LC_ALL=C \
/home/yanghj/my_data/software/container/TETools202309.sif BuildDatabase -name MHP25F ../MHP25F_Chr.fa
# 运行 RepeatModeler Repeatmodeler是通过重续序列的结构特征来进行从头注释
singularity exec --env SINGULARITYENV_LANG=C --env SINGULARITYENV_LC_ALL=C \
/home/yanghj/my_data/software/container/TETools202309.sif RepeatModeler \
-database MHP25F \
-threads 50 -LTRStruct 1>repeatmodeler.log 2>&1
# 运行 RepeatMasker Repeatmasker 基于与已知的重复序列数据库比对来寻找重复序列，
singularity exec --env SINGULARITYENV_LANG=C --env SINGULARITYENV_LC_ALL=C \
/home/yanghj/my_data/software/container/TETools202309.sif RepeatMasker \
-e rmblast -pa 50 -qq \
-lib MHP25F-families.fa ../MHP25F_Chr.fa 1>repeatmasker.log 2>&1
# RepeatMasker自带程序 将out文件转化成bed格式
# conda acivate EDTA环境 最后转化后的文件MHP25F_Chr.fa_rm.bed
RM2Bed.py genome.fa.out
# MHP25F_Chr.fa.out 注释gff文件，转化成bed格式
# MHP25F_Chr.fa.tbl 注释统计表格
TETools202309.sif文件
通过百度网盘分享的文件：TETools202309.sif
链接：https://pan.baidu.com/s/1QmKURNpjmLPRE0HZTre_5g
提取码：8pk8
方法三：用EDTA进行重复序列注释，会花费时间很久
gtihub网址：https://github.com/oushujun/EDTA
conda create -n EDTA
conda activate EDTA
mamba install -c conda-forge -c bioconda edta
# EDTA安装的 repeatMasker不能直接使用，缺一个数据库
# 最好提供了RepBase数据，百度云盘下载
通过百度网盘分享的文件：RepBase29.03.fasta.tar.gz
链接：https://pan.baidu.com/s/1yoHTlJUYZ-t4ZLmrk64U8Q
提取码：otot
tar -zxvf RepBase29.03.fasta.tar.gz
tar -zxvf RepBaseRepeatMaskerEdition-20181026.tar.gz
添加数据库
(EDTA) [yanghj@cln01 RepeatMasker]$ ./addRepBase.pl -libdir Libraries/
Rebuilding RepeatMaskerLib.h5 master library
- Read in 49011 sequences from Libraries//RMRBSeqs.embl
- Read in 49011 annotations from Libraries//RMRBMeta.embl
Merging Dfam + RepBase into RepeatMaskerLib.h5 library...
Database: Dfam withRBRM
Version: 3.2
Date: 2020-07-02
Dfam - A database of transposable element (TE) sequence alignments and HMMs.
RBRM - RepBase RepeatMasker Edition - version 20181026
Total consensus sequences: 51780
Total HMMs: 6915
大家在这里，请不要直接复制，结合自己的文件进行调整。
# 进行重复序列注释的代码
EDTA.pl --genome genome.fa --species others --overwrite 0 --sensitive 1 --evaluate 1 --anno 1 --threads 100
# --sensitive 默认是0，修改为1，调用repeatmodeler，花时间更长， 是否用RepeatModeler分析剩下的TE：0
# --evaluate 启用评估模式，将评估结果输出到文件中
# --anno 指定进行转座子注释，"1"代表进行完整的注释
若我们的教程对你有所帮助，请点赞+收藏+转发，大家的支持是我们更新的动力！！
往期部分文章
1. 最全WGCNA教程（替换数据即可出全部结果与图形）
推荐大家购买最新的教程，若是已经购买以前WGNCA教程的同学，可以在对应教程留言，即可获得最新的教程。（注：此教程也仅基于自己理解，不仅局限于此，难免有不恰当地方，请结合自己需求，进行改动。）
2. 精美图形绘制教程
3. 转录组分析教程
4. 转录组下游分析
小杜的生信筆記 ，主要发表或收录生物信息学教程，以及基于R分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料!!
