参考0
参考1
miRNA预测工具miRDeep-P2简明教程
植物miRNA靶基因预测
下机数据整理汇总
find -name ./msj *.clean.fastq.gz |xargs -i cp {} ./msj/clean/
查找所有的质控过的数据，移动到clean文件夹。
我的是水稻的miRNA数据。所以先下载水稻的各种文件。
下载地址
包括基因组序列、基因组注释、基因组蛋白质注释、基因组cds序列。
下载miRNA数据
下载miRNA注释文件
1. 测序下机数据质控、去接头、检测分布
1.1 质控fastqc
fastqc -i ./clean/*.fastq.gz -t 4 #4个进程同时运行
1.2 去接头
miRNA的一般用cutadapt,同时去掉reads中的adapter，低质量的reads以及过长过短的reads。
cutadapt下载安装配置环境变量，参考官网
环境变量在~/.bashsrc_profile
去掉接头。
cutadapt -a AGATCGGAAGAG --quality-base 33 -m 10 -q 20 --discard-untrimmed -o $Co1.fa >cutadapt.info
具体的接头文件，在fastqc的输出文件中，找adapter可以看到。
1.3 检测片段分布
新建length-distribute.py
python 代码如下：
#!/usr/bin/env/ python3
#Author:Frank
import sys
miRNA_len={}
#此处的60长度根据你的每行数据的最长长度设置，过小可能会报错。
for i in range(0,60):
miRNA_len[i]=0
for i in open(sys.argv[1]):
if i.startswith('@') or i.startswith('+'):
continue
length=len(i)-1
miRNA_len[length] += 1
for i in miRNA_len:
print(str(i)+"\t"+str(miRNA_len[i]/2))
使用方法：
示例：python3 length-distribute.py ./merge/S168-Na-1_L1.fastq >S168-1.stat
实际使用脚本循环所有数据：
2. 拼接
miRNA拼接最好使用bowtie，而不是bowtie2.详情地址
2.1 bowtie的使用
-建立索引
bowtie-build -f genome.fasta genome.fa &
bowtie-build -f miRNA.fa miRNA.fa &
bowtie -q -v 2 -l 10 -k 15 genome.fa C01.fa -S c01.sam 2>mapping.info
此处比对我实际使用的脚本循环，因为样本量太多了。
name_array=(Co1 S168-Na STTM168 Co1-Na)
for ((n=0;n<4;n++));
do
for ((i=1;i<4;i++));
do
read="./merge/"${name_array[$n]}"-"$i"_L1.fastq"
bowtie -q -v 2 -l 10 -k 15 genome.fa $read -S ${name_array[$n]}"-"$i".sam" 2>>mapping.info
done
done
在输出文件mapping.info中可以找到mapping到的比率。
整体在90%以上，数据比较好。
对于miRNA比对一般有2种方式，一种是比对到基因组上，另一种是比对到对应物种的miRNA数据库中。此处的genome.fa是基因组文件。
-mi.fa 是把miRNA.fa的U碱基转换成T碱基。
-miRNA.fa 是原始的从mirbase下载的茎环结构的数据，自己命名。
bowtie分别构建上述三个基因组的索引文件。
参考曾建明多年前的 提问，miRNA比对到mirBase的miRNA数据库时，是需要把U转换为T的。
但是事实上，转变碱基之后的比对效果仍然很低（0.1-0.2%）。不知道原因。所以我暂时使用的是和基因组的比对结果。
3. 安装使用HTSeq（非root用户）（据说可以使用featurecount替代HTSeq）
作为非root用户，安装许多软件非常头疼。总是没有权限。
先在GitHub下载htseq。
下载地址
unzip htseq-master
python setup.py build install --user
安装完成，添加环境变量。
vim ~/.bash_profile
export PATH=$HOME/chaim/.local/bin:$PATH
重启终端，测试htseq。
htseq-count -h
会出现相应的帮助信息。
HTSeq-count使用API
htseq统计counts的脚本
#!/bin/bash
name_array=(Co1 S168-Na STTM168 Co1-Na)
for ((n=0;n<4;n++));
do
for ((num=1;num<4;num++));
do
htseq-count -s no -t miRNA -i ID -o ${name_array[$n]}"-"$num".hc.sam" ${name_array[$n]}"-"$num".sam" miRNA.gff3 | tee ${name_array[$n]}"-"$num".count"
ls ${name_array[$n]}"-"$num".count" >> count.list
#echo ${name_array[$n]}"-"$num"count"
done
done
吐槽一下,水稻基因组本身是有2种命名格式的，基因注释是有Chr1和Chr01的格式的。所以一定要看清楚自己的注释文件和Sam文件的格式是否一致，否则会出现counts=0.
4. 表达量差异分析（R语言中操作）
无生物学重复使用DEGseq，有生物学重复使用DESeq2
无生物学重复DEGseq（以下代码摘自引用1，未经过实践。）
library("DEGseq")
#BT_0_1
geneExpMatrix1<- readGeneExp("ht.genotype_data1.txt", geneCol=1, valCol=3)
geneExpMatrix2<- readGeneExp("ht.genotype_data2.txt", geneCol=1, valCol=2)
write.table(geneExpMatrix1[30:31,],row.names=FALSE)
write.table(geneExpMatrix2[30:31,],row.names=FALSE)
pdf(file="data1_2.pdf")
layout(matrix(c(1,2,3,4,5,6),3, 2, byrow=TRUE))
par(mar=c(2,2, 2, 2))
DEGexp(geneExpMatrix1=geneExpMatrix1,geneCol1=1, expCol1=2, groupLabel1="data1",geneExpMatrix2=geneExpMatrix2,geneCol2=1, expCol2=2,groupLabel2="data2",method="MARS",outputDir="05DEmiRNA/DEGSeq")
dev.off()
有生物学重复DESeq2 (RNA-seq也可以使用它计算表达量)
网友评论