[【基因組注釋】ncRNA注釋 原創(chuàng)]（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
手動(dòng)解壓

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

將兩個(gè)數(shù)據(jù)庫得到的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標(biāo)準(zhǔn)序列

參考http://www.reibang.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}' 輸出每個(gè)項(xiàng)目第五個(gè)的內(nèi)容
uniq 命令刪除重復(fù)的行
ensemble 文件中的RNA類型

圖片10.png

將其中某種類型的RNA（Mt_tRNA）提出來放入單個(gè)文件中

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

1.從sra文庫中下載fastq 文件

2.fastqc 檢查fastq質(zhì)量文件

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

去除低質(zhì)量堿基，使用作者給的-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

查看以上每個(gè)文件的行數(shù), 就可以繪制下表鸡号。每4行代表一個(gè)reads

image.png

 wc -l DRR34.clean.fastq

image.png

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

利用其他人寫的pl文件

perl groupReads.pl DRR34.clean.fastq

生成以下三個(gè)文件

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 作為輸入文件，大大提高運(yùn)行速度

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個(gè)文件临梗，一個(gè)是匹配了的fasta文件，一個(gè)是未匹配成功的fasta文件序列。還一個(gè)匹配生的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扣典，所有得替換负甸，只能用正則表達(dá)是流强。目前還沒想到
但是前面匹配過程中已經(jīng)給出了匹配得fasta 序列，如 DRR34.genome.fasta

6 MIRbase 匹配

https://mirbase.org/ 進(jìn)入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個(gè)文件

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

構(gòu)建索引

makeblastdb -in hsa.Rfam.fasta -dbtype nucl

image.png

進(jìn)行比對

 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

該結(jié)果與作者展示的結(jié)果不一致，是由于Rfam的版本不一樣蚕捉，我選擇的是14.8版本奏篙，而作者選用的是11的版本

image.png

11版本的結(jié)果

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注釋結(jié)果

image.png

改為11版本的注釋就行

8 不同miRNA[diffential expression]

根據(jù)之前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

四個(gè)矩陣進(jìn)行合并

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

image.png