[【基因组注释】ncRNA注释 原创]（https://blog.51cto.com/u_15668923/5639934）
http://www.chenlianfu.com/?p=2901

image.png

image.png

image.png

image.png

image.png

1. RFAM.fa 文件

从rfam下载所有的序列##

wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.8/fasta_files/*
cat RF*.gz>Rfam.all.fa.gz
手动解压

Convert a multi-line FASTA file to a single-line FASTA file:

cat Rfam.all.fa | fasta_formatter -w 0 >Rfam.nowrap.fa

提取fasta 文件中的rRNA 序列

grep -A 1 "rRNA" Rfam.nowrap.fa > hg19.rRNA.fa

将两个数据库得到的fa合并在一起

cat hg19.rRNA.fa rRNA.${species}.fa > all.human.rRNA.fa

2. 从ensemble下载ncRNA.fa

Convert a multi-line FASTA file to a single-line FASTA file:

cat Homo_sapiens.GRCh38.ncrna.fa| fasta_formatter -w 0 > Rfam.nowrap.fa

3 下载mt-tRNA标准序列

参考https://www.jianshu.com/p/ef20c4932b3b

$ wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/ncrna/Homo_sapiens.GRCh38.ncrna.fa.gz
$ gunzip Homo_sapiens.GRCh38.ncrna.fa.gz 
$ head -100 Homo_sapiens.GRCh38.ncrna.fa #查看fa文件前100行

查看fa文件小RNA类型

cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/^>/p' | awk '{print $5}' | uniq

ENST00000459236.1 ncrna chromosome:GRCh38:7:128443449:128443510:-1 gene:ENSG00000238590.1 gene_biotype:snRNA transcript_biotype:snRNA gene_symbol:RNU7-54P description:RNA, U7 small nuclear 54 pseudogene [Source:HGNC Symbol;Acc:HGNC:34150]
TAGTGTTACAGCGCTTTGAGAATTTGTCTAGCAGGCCTTTGGGTCTTTACTGGAAAACCC
CT

sed 找到所有以>开头的行
awk '{print $5}' 输出每个项目第五个的内容
uniq 命令删除重复的行
ensemble 文件中的RNA类型

图片10.png

将其中某种类型的RNA（Mt_tRNA）提出来放入单个文件中

mkdir rfam
cd rfam
touch Mt_tRNA.fa ##新建文件  
cd  在返回上一级
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/Mt_tRNA/{p;n;p}' >>  ./rfam/Mt_tRNA.fa

1.从sra文库中下载fastq 文件

2.fastqc 检查fastq质量文件

3. cutadapt 去除

#去除有5adaptor污染的reads
cutadapt -g GTTCAGAGTTCTACAGTCCGACGATC --discard-trimmed DRR34.fastq -o DRR34.5adptor.fastq
#去除3adaptor
cutadapt -a TCGTATGCCGTCTTCTGCTTG -m 17 -M 35 --overlap=5 --discard-untrimmed DRR34.5adptor.fastq -o DRR34.3adptor.fastq

去除N的碱基

./removeN DRR34.3adptor.fastq DRR34.unknowNul.fastq

去除低质量碱基，使用作者给的-Q 33 不行，改为64可以

fastq_quality_filter -q 20 -p 80 -i DRR34.unknowNul.fastq -o DRR34.clean.fastq  -Q 64
~/Biosofts/fastqc/FastQC/fastqc DRR34.clean.fastq -t 2 -o qcOutdir

image.png

查看以上每个文件的行数, 就可以绘制下表。每4行代表一个reads

image.png

 wc -l DRR34.clean.fastq

image.png

4.[查看长度分布length distribution]

利用其他人写的pl文件

perl groupReads.pl DRR34.clean.fastq

生成以下三个文件

image.png

image.png

image.png

image.png

利用写好的R程序生成图

Rscript barplot.R DRR34.lengthDistribution.txt DRR34.lengthDistribution.pdf

image.png

5 bowtie 基因组mapping

使用上述的uniq.fasta 作为输入文件，大大提高运行速度

bowtie -f -p 4 -x ./human_index DRR34.uniq.fasta -S DRR34.clean.sam --un DRR34.clean.
unaligned.ncRNA.fa --al DRR34.genome.fasta

生成3个文件，一个是匹配了的fasta文件，一个是未匹配成功的fasta文件序列。还一个匹配生的SAM文件

image.png

samtools view -hSb  DRR34.clean.sam -o DRR34.clean.bam
samtools view -bS DRR34.clean.sam | samtools sort -o DRR34.clean.bam 
 samtools flagstat DRR34.clean.bam

image.png

查看sam文件格式

less -SN DRR34.clean.sam

image.png

image.png

染色体不为*则代表匹配成功的序列
需要筛选除匹配上的序列

awk '$3~/chr/{print ">"$1"\n"$1}' DRR34.sam >DRR34.genome.fasta

但是这里的染色体开头并没有chr，所有得替换，只能用正则表达是。目前还没想到
但是前面匹配过程中已经给出了匹配得fasta 序列，如 DRR34.genome.fasta

6 MIRbase 匹配

https://mirbase.org/ 进入download wget下载

image.png

image.png

perl pickSpeciesRNA.pl --input mature.fa --out hsa.mature.fa --species 'Homo sapiens'
perl pickSpeciesRNA.pl --input hairpin.fa --out hsa.hairpin.fa --species 'Homo sapiens'

建立blast 索引库

makeblastdb -in hsa.mature.fa -dbtype nucl
makeblastdb -in hsa.hairpin.fa -dbtype nucl

开始blast

blastn -query DRR34.genome.fasta -db hsa.mature.fa -task blastn-short -outfmt 6
 -evalue 0.01 -out DRR34.mature.blast.txt -num_alignments 1

生成得DRR34.mature.blast.txt文件格式如下

image.png

同样
blastn -query DRR34.genome.fasta -db hsa.hairpin.fa -task blastn-short -outfmt 6 -evalue 0.01 -out DRR34.hairpin.blast.txt -num_alignments 1

perl miRNAStat.pl --mature DRR34.mature.blast.txt --hairpin DRR34.hairpin.blast.txt --genomeFa DRR34.genome.fasta 
 --readCount DRR34.readCounts.txt --fastaOut DRR34.removeMirna.fasta --statOut DRR34.ncRNA.stat.txt --expOut DRR34.exp.txt --mismatch 1

生成3个文件

image.png

image.png

image.png

image.png

7 其他coding RNA 匹配

image.png

perl pickSpeciesRNA.pl --input Rfam.fa --out hsa.Rfam.fasta --species 'Homo sapiens'

image.png

构建索引

makeblastdb -in hsa.Rfam.fasta -dbtype nucl

image.png

进行比对

 blastn -query DRR34.removeMirna.fasta -db hsa.Rfam.fasta -task blastn-short 
-outfmt 6 -evalue 0.01 -out DRR34.Rfam1.blast.txt -num_alignments 1

image.png

image.png

该结果与作者展示的结果不一致，是由于Rfam的版本不一样，我选择的是14.8版本，而作者选用的是11的版本

image.png

11版本的结果

image.png

perl RfamStat.pl --Rfam DRR34.Rfam1.blast.txt --inputFa DRR34.removeMirna.fasta --readCount DRR34.readCounts.txt --fastaOut DRR34.removeRfam.fasta --statOut DRR34.ncRNA.stat.txt

image.png

最终各类RNA注释结果

image.png

改为11版本的注释就行

8 不同miRNA[diffential expression]

根据之前perl miRNAStat.pl --mature DRR34.mature.blast.txt --hairpin DRR34.hairpin.blast.txt --genomeFa DRR34.genome.fasta
--readCount DRR34.readCounts.txt --fastaOut DRR34.removeMirna.fasta --statOut DRR34.ncRNA.stat.txt --expOut DRR34.exp.txt --mismatch 1
得到的DRR34.exp.txt文件，同理其他文件也是类似得到

image.png

四个矩阵进行合并

perl getExpMatrix.pl DRR33.exp.txt DRR34.exp.txt DRR35.exp.txt DRR36.exp.txt

image.png