miRNA分析--比對（二）

 miRNA分析--數據過濾（一）

 

在比對之前為了減少比對時間，將每一個樣本中的reads進行合並，得到fasta格式，其命名規則如下：

樣本_r數子_x數字

r 中的數字表示reads序號；

x 中的數字表示該條reads重復次數

 

比對分為兩條策略

1、根據本物種已有的miRNA序列進行比對，

已知當miRNA序列從 miRBase 或者 sRNAanno得到

 

(應該將clean reads比對到所研究物種到tRNA, rRNA, snoRNA,mRNA等數據，允許一個錯配，將比對上等reads過濾，也可以比對到參考基因組，將為未比對到到reads過濾掉，但是本次我沒有這么做)

對於第一種情況，我采用bowtie將reads比對到成熟miRNA

##建立索引
bowtie-build ref.fa
##比對
bowtie -v 0 -m 30 -p 10 -f  ref.fa sample.fa sample.bwt
參數解釋
-v: 允許0個錯配
-p: 10個線程 
-m: 當比對超過這個數時，認為時未比對
-f： 輸入序列fasta

根據.bwt 文件可以計算出每個已知當miRNA中比對上當reads數量，別忘記乘以 x后面的數

2、直接比對到參考基因組並進行新miRNA鑒定

采用miR-PREFeR進行Novel miRNA鑒定

其githup主頁：https://github.com/hangelwen/miR-PREFeR

• 需安裝ViennaRNA ( 最好是1.8.5或2.1.2, 2.1.5版本)

#我裝的是最新版
wget https://www.tbi.univie.ac.at/RNA/download/sourcecode/2_4_x/ViennaRNA-2.4.10.tar.gz 
tar zvxf ViennaRNA-2.4.10.tar.gz
cd ViennaRNA-2.4.10
./configure --prefix="/user/tools/ViennaRNA/" --without-perl
make
make install

• 安裝

 git clone https://github.com/hangelwen/miR-PREFeR.git 

• 數據准備

1⃣️ref.fa

2⃣️miRNA 比對到ref.fa的sam文件 （sam 文件中的reads 必須是collapse reads ）

3⃣️gff 文件（可選，記錄需要屏蔽掉的信息，比如重復序列等）

bowtie 比對

 bowtie -a -v 0 -m 30 -p 10 -f ref.fa sample.fa -S sample.sam 

准備configure 文件

#example configuration file for the miR-PREFeR pipeline.
#lines start with '#' are comments

#miR-PREFeR path, please change to your path to the script folder. 
#Absolute path perfered.
PIPELINE_PATH = /miR

#Genomic sequence file in fasta format.  Absolute path perfered. If a path
#relative if used, it's relatvie to the working directory where you execute
#the pipeline.
FASTA_FILE =  genome_v1.fa

#Small RNA read alignment file in SAM format. The SAM file should contain 
#the SAM header. If N samples are used, then N file names are listed here, 
#separated by comma. please note that before doing alignment, process the 
#reads fasta files using the provided script 'process-reads-fasta.py' to
#collapse and rename the reads. Absolute path perfered. If a path
#relative if used, it's relatvie to the working directory where you execute
#the pipeline.
ALIGNMENT_FILE = ./trm_XX-1_L1_I309.R1.fastq_trm_fa.fa.sam, ./trm_XX-2_L1_I310.R1.fastq_trm_fa.fa.sam, ./trm_XX-3_L1_I311.R1.fastq_trm_fa.fa.sam, ./trm_XY-1_L1_I312.R1.fastq_trm_fa.fa.sam, ./trm_XY-2_L1_I313.R1.fastq_trm_fa.fa.sam, ./trm_XY-3_L1_I314.R1.fastq_trm_fa.fa.sam, ./trm_YY-1_L1_I315.R1.fastq_trm_fa.fa.sam, ./trm_YY-2_L1_I316.R1.fastq_trm_fa.fa.sam, ./trm_YY-3_L1_I332.R1.fastq_trm_fa.fa.sam

#GFF file which list all existing annotations on genomic sequences FASTA_FILE.
#If no GFF file is availble, comment this line out or leave the value blank.
#Absolute path perfered. If a path relative if used, it's relatvie to the
#working directory where you execute the pipeline.
#CAUTION: please only list the CDS regions, not the entire miRNA region, because 
#miRNAs could be in introns. This option is mutual exclusive with 'GFF_FILE_INCLUDE'
#option.
# If you have a GFF file that contains regions in which you want to predict whehter 
# they include miRNAs, please use the 'GFF_FILE_INCLUDE' option instead.
GFF_FILE_EXCLUDE = CDS.gff

# Only predict miRNAs from the regions given in the GFF file. This option is mutual 
# exclusive with 'GFF_FILE_EXCLUDE'. Thus, only one of them can be used.
#GFF_FILE_INCLUDE = ./TAIR10.chr1.candidate.gff

#The max length of a miRNA precursor. The range is from 60 to 3000. The default
#is 300. 
PRECURSOR_LEN = 300

#The first step of the pipeline is to identify candidate regions of the miRNA
#loci. If READS_DEPTH_CUTOFF = N, then bases that the mapped depth is smaller
#than N is not considered. The value should >=2.
READS_DEPTH_CUTOFF = 20

#Number of processes for this computation. Using more processes speeds up the computation,
#especially if you have a multi-core processor. If you have N cores avalible for the 
#computation, it's better to set this value in the range of N to 2*N.
#If comment out or leave blank, 1 is used. 
NUM_OF_CORE = 4

#Outputfolder. If not specified, use the current working directory. Please make sure that 
#you have enough disk space for the folder, otherwise the pipeline may fail.
OUTFOLDER = spinach-result

#Absolute path of the folder that contains intermidate/temperary files during the
# run of the pipeline. If not specified, miR-PREFeR uses a folder with suffix "_tmp"
#under OUTFOLDER by default. Please make sure that you have enough disk space for the
# folder, otherwise the pipeline may fail.
#TMPFOLDER = /tmp/exmaple
TMPFOLDER = 

#prefix for naming the output files. For portability, please DO NOT contain any 
#spaces and special characters. The prefered includes 'A-Z', 'a-z', '0-9', and 
#underscore '_'.
NAME_PREFIX = spinach-example

#Maximum gap length between two contigs to form a candidate region. 
MAX_GAP = 100

# Minimum and maximum length of the mature sequence. Default values are 18 and 24.
MIN_MATURE_LEN = 18
MAX_MATURE_LEN = 24

# If this is 'Y', then the criteria that requries the star sequence must be expressed 
# is loosed if the expression pattern is good enough (.e.g. the majority of the reads 
# mapped to the region are mapped to the mature position.). There are lots of miRNAs 
# which do not have star sequence expression. The default value is Y.
ALLOW_NO_STAR_EXPRESSION = Y

# In most cases, the mature star duplex has 2nt 3' overhangs. If this is set to 'Y', then
# 3nt overhangs are allowed. Default is 'N'.
ALLOW_3NT_OVERHANG = N

#The pipeline makes a checkpoint after each major step. In addition, because the
#folding stage is the most time consuming stage, it makes a checkpiont for each
#folding process after folding every CHECKPOINT_SIZE sequences. If the pipeline
#is killed for some reason in the middle of folding, it can be restarted using
#'recover' command from where it was stopped. The default value is 3000. On my
#system this means making a checkpoint about every 5 minutes.
CHECKPOINT_SIZE = 3000

運行

 python miR_PREFeR.py -L -k pipeline configfile 

pipeline里包含prepare, candidate, fold, predict四步。如果某步中斷了，還可以續跑

 python miR_PREFeR.py -L recover configfile 

 

輸出結果

根據example_miRNA.detail.csv 文件 寫一個腳本 提取每個miRNA的reads 數量，進而做差異分析

 

差異分析

差異分析采用DESeq2， 可看我之前寫的miRAN 分析以及mRNA分析

 

 

 

關注下方公眾號可獲得更多精彩

 

ref

1、計算已知miRNA當表達量

2、省心省事的植物miRNA分析軟件miR-PREFeR,值得擁有