miRNA数据分析流程

一般流程：参考miRNA-seq分析流程 - 简书

一篇文章学会miRNA-seq分析 - 云+社区 - 腾讯云

质控:FastQC

去接头：trim_Galore去除Illumina接头；cutadapt去除自定义接头

参考基因比对：两种，一种是和参考基因组比对，可以得到新的未预测的reads；另一种是和已知数据库比对，操作简单。.对于50bp以下的样本用bowtie，50-200bp的样本用bowtie2，二者区别参考以下链接bowtie和bowtie2使用条件区别及用法_soyabean555999的博客-CSDN博客_bowtie2

定量：Featurecounts/ HTseq。featurecounts的使用说明 - 简书

我们本次分析是先进行bedtools interact找参考基因和mapping reads的交集，后用python计数。

表达量差异分析（R语言中做）：无生物学重复使用DEGseq，有生物学重复使用DESeq2

miRNA样本配对mRNA表达量获取

mi下游数据分析表达

​

本次数据分析流程： 

bedtools interact 参考网址：bedtools intersect用法 (intersectBed) - 简书

                                             bedtools intersect 的八个常用案例 - 简书 (jianshu.com)

定量：FeatureCounts，但出问题了，匹配率过低，只有0.3%左右。

           原因：1、参考GTF的问题

                     2、rRNA contimination

featureCount的速度很快，但输出文件中只有数据统计，没有reads注释信息，如果再进行注释，需要进行bedtools annotation或者使用python根据gene position提取GTF中的信息。但GTF也很大。

目前存在的最大问题：三种计数方法得出的结论不一致。

1、bowtie后得到sam文件，转为bam文件后，进一步转为bed文件，直接用python计数（用之前tRNA的code即可）

2、bowtie后进行bedtools interact 寻找交集，输出A文件与B文件相交的信息，相当于得到了一个未计数但注释了的文件。进一步使用python计数即可，但该文件非常大，运行需要大概6h得到一份数据，注意提取I列时需要先统计有多少行。运行完后得到了map的gene的position及ID，counts等信息，进一步使用excel的VLOOKUP函数提取gene description，gene description的文件下载参考此文章参考基因组和基因注释文件 | Public Library of Bioinformatics，使用的下载网站为BioMart http://www.ensembl.org/biomart/martview。 至此得到了注有基因功能及name，counts的表格，但由于不同网站对基因的命名及ID编号信息不一致，导致提取到的gene description信息不完全。如果想要解决此问题，需要重新用相同网站ensembl下载的GTF进行bedtools interact。最后为了得到数据库之外的新型RNA，需与miRNA database里已有的数据进行比较，以sequence为参考进行比较，则使用bedtools getfasta进行sequence提取，提取方法：将chr,gene start, gene stop三列从excel中复制到记事本中，在linux中直接将文件后缀改为.bed文件，进一步使用getfasta即可。

得到后写code用python比较即可。

3、FeatureCounts，但出问题了，匹配率过低，只有0.3%左右。

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章目录

• 前言
• 一、pandas是什么？
• 二、使用步骤
  • 1.引入库
  • 2.读入数据
• 总结

---

前言

       随着测序技术的发展，数据分析这门技术也越来越重要，作为生信小白自学着实有些艰难，分享一篇学习miRNA数据分析的记录，数据分析结果是否可靠仍然有待考究。

---

一、miRNA是什么？

MicroRNAs (miRNAs) are small RNA molecules, which are ∼22 nt sequences that have an important role in the translational regulation and degradation of mRNA by the base's pairing to the 3′-untranslated regions (3′-UTR) of the mRNAs. The miRNAs are derived from the precursor transcripts of ∼70–120 nt sequences, which fold to form as stem–loop structures, which are thought to be highly conserved in the evolution of genomes. Previous analyses have suggested that ∼1% of all human genes are miRNA genes, which regulate the production of protein for 10% or more of all human coding genes.

二、使用步骤



1.数据下载及质控

代码如下：

$ conda activate rnaseq #打开软件安装的环境，实验室电脑安装在这个环境之下，所以每次使用前先打开
$ fastqc #该软件可视化，打开后就可以看到了

2.trim_galore 

代码如下：

$ trim_galore -q=20 --phred33 --fastqc --illumina --length 20 --stringency 1 -e 0.1 -o file inputfile

$cutadapt -a CCCAGATCGGAAGAG -g AAAAAAAAAAAAAAAAAAAA --error-rate=0.0 --times=2 --overlap=3 --minimum-length=20 --output=/home/caolab/PT/Data_analysis/clean_data_after_cutadapt/caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_cutadapted.fq.gz /home/caolab/PT/Data_analysis/cleandata_after_trim_galore/caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_trimmed./caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_trimmed.fq.gz

 3.cutadapt 

代码如下：cutadapt 使用详解 - 简书

Usage:
    cutadapt -a ADAPTER options input.fastq
    cutadapt -a AACCGGTT input.fastq > output.fastq
For paired-end reads:
    cutadapt -a ADAPT1 -A ADAPT2 [options] -o out1.fastq -p out2.fastq in1.fastq in2.fastq
$ cutadapt -a CCCAGATCGGAAGAG -g AAAAAAAAAAAAAAAAAAAA --error-rate=0.0 --times=2 --overlap=3 --minimum-length=20 --output=/home/caolab/PT/Data_analysis/clean_data_after_cutadapt/caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_cutadapted.fq.gz /home/caolab/PT/Data_analysis/cleandata_after_trim_galore/caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_trimmed./caobo20210822-S10_FKDL210228552-1a-D710-AK1546_2_trimmed.fq.gz
常用参数
-g: #剪切reads 5'端adapter(双端测序第一条read)，加$表示adapter锚定在reads 5'端
-a: #剪切reads 3'端adapter(双端测序第一条read)，加$表示adapter锚定在reads3'端
-O MINLENGTH, --overlap=MINLENGTH #adapter与reads最小overlap,才算成功识别; Default: 3
-m LENGTH, --minimum-length=LENGTH: 根据最短长度筛选reads;Default: 0
--discard-untrimmed, --trimmed-only #丢掉不包含adapter的reads
 -e ERROR_RATE, --error-rate=ERROR_RATE  #adapter匹配允许的最大错配率（错配/匹配片段长度)；Default: 0.1
--no-trim: 不修剪adapter，直接输出满足跳进啊的reads
-u LENGTH, --cut=LENGTH:  #修剪reads 5'/3'端碱基,正数：从开始除移除碱基；负数：从末尾处移除碱基；
-q [5'CUTOFF,]3'CUTOFF, --quality-cutoff=[5'CUTOFF,]3'CUTOFF: #修剪低质量碱基
-l LENGTH, --length=LENGTH: #将reads修剪的最终长度
--trim-n: #修剪reads末端的'N'
-o FILE, --output=FILE: #输出文件
--info-file=FILE：每条reads和与reads匹配的adapter的信息
--too-short-output=FILE: #为reads长度最小值设定阈值筛选reads后，要丢弃的部分输出到文件；长度依据m值设定；   
--too-long-output=FILE：#为reads长度最大值设定阈值筛选reads后，要丢弃的部分输出到文件；长度依据M值设定； 
--untrimmed-output=FILE: #将没有adapter未做修剪的reads输出到一个文件;默认输出到trimmed reads结果文件
--max-n=COUNT：#reads中N的数量，设定整数或小数(N的占比)

双端测序参数
-A ADAPTER：  #第二条reads 3'adapter
-G ADAPTER：#第二条reads 5'adapter
-U LENGTH： #从第二条reads上修剪的长度
-p FILE, --paired-output=FILE： #第二条reads的输出结果
--untrimmed-paired-output=FILE：#第一条reads没有adapter时，将第二条reads输出到文件；默认输出到trimmed reads结果文件   

 4.bowtie2 

代码如下：

usage：
$ bowtie2 [options]* -x  {-1  -2  | -U } [-S 
实例
$ bowtie2 -x /home/caolab/PT/Data_analysis/S1-5938.989kb.seq.fasta -U /home/caolab/PT/Data_analysis/clean_data_after_cutadapt/caobo20210822-S2_FKDL210228552-1a-D702-AK1543_2_cutadapted -S /home/caolab/PT/Data_analysis/sam_after_bowtie2/caobo20210822-S2_FKDL210228552-1a-D702-AK1543_2_bowtie2
#bowtie2之前需要给参考基因组建立索引

 5. samtools sort

代码如下：

usage：
$ samtools sort caobo20210822-S2_FKDL210228552-1a-D702-AK1543_2_bowtie2 > caobo20210822-S2_FKDL210228552-1a-D702-AK1543_2_bowtie2.bam
#将bowtie2 输出的bam文件转为bam文件,同时进行排序

  6.bedtools bamtobed 

代码如下：

bamtobed
bedtools bamtobed -i /run/media/caolab/B47415C574158AEE/E盘/microRNA/上海分析/miRNA_result/bam/GT_26_9_M_DKDL202002780-1a_H7YJMCCX2_L3_sam_sorted.bam > /home/caolab/miRNA/bed/GT_26_9_M_DKDL202002780-1a_H7YJMCCX2_L3_sam_sorted.bed
#bedtools中用>表示输出文件

  7.bedtools getfasta

代码如下：

$ bedtools getfasta -fi /home/caolab/miRNA/hg38.fa -bed /home/caolab/miRNA/bed/GT_26_9_M_DKDL202002780-1a_H7YJMCCX2_L3_sam_sorted.bed -fo /home/caolab/miRNA/extract_DNA_sequence/ GT_26_9_M_DKDL202002780-1a_H7YJMCCX2_L3_sam_sorted.fa
#提取sequence，对应galaxy的extract genomic DNA

  8. featureCounts

代码如下：featurecounts的使用说明 - 简书

install featureCounts:
$ conda activate rnaseq
$ wget https://jaist.dl.sourceforge.net/project/subread/subread-1.6.0/subread-1.6.0-Linux-x86_64.tar.gz#下载软件的命令
$ tar -zxvf subread-1.6.0-Linux-x86_64.tar.gz#对下载的软件进行解压的命令，注意不要打错字母，否则会一直报错。
$ subread-1.6.0-Linux-x86_64/bin/featureCounts#打开featureCounts的命令

usage：
$ featureCounts [options] -a  -o  input_file1 [input_file2]
实例：
$ subread-1.6.0-Linux-x86_64/bin/featureCounts -T 10 -a /home/caolab/miRNA/gencode.v29.annotation.gtf/gencode.v29.annotation.gtf -t exon -g gene_id -o counts.txt /home/caolab/miRNA/bam/GT_26_9_M_DKDL202002780-1a_H7YJMCCX2_L3_sam_sorted.bam

 9. 用python对bedtools interact 之后的文件进行计数，并提取DEI列，提取DE两列为gene position，I列为name，ID等

代码如下：thinkbook笔记本上的python在E盘中，相关code也在里面。做这个时要用打电脑进行，笔记本带不起来，大电脑的windows中的pycharm在D盘的E盘文件夹中。

code需要从大电脑上copy过来

10. 用python对bed文件进行计数，即下文中的###统计数量###部分（先溢）

代码如下：

import xlrd
import xlwt
from xlutils.copy import copy
date={}


#####################################统计数量###############################

def xlsxdate(date):
    writebook = xlwt.Workbook()  # 打开一个excel
    sheet = writebook.add_sheet('test')  # 在打开的excel中添加一个sheet
    i = 0
    sheet.write(i, 0, '名字')  # 写入名字
    sheet.write(i, 1, '开始位置')  # 开始位置
    sheet.write(i, 2, '结束位置')  # 结束位置
    sheet.write(i, 3, '数量')  # 数量
    sheet.write(i, 4, '正负')  # +=
    for s in date:
        if date[s]['sum'] >= 10:
            i += 1
            sheet.write(i, 0, date[s]['name'])  # 写入excel，i行0列
            sheet.write(i, 1, date[s]['start'])
            sheet.write(i, 2, date[s]['end'])  # 写入excel，i行0列
            sheet.write(i, 3, date[s]['sum'])
            sheet.write(i, 4, date[s]['e'])  # 写入excel，i行0列
    writebook.save('date2.xlsx')  # 一定要记得保存
    print(i)
    pass

with open('Galaxy398-[Filter_on_data_319__Filtered_BAM_(as_BED)].bed') as f:

    for line in f:
       name,start,end,c,d,e=line.split()
       s=name+start+end
       if s in date:
           date[s]['sum']+= 1
       else:
           date[s]={}
           date[s]['name']=name
           date[s]['start']=start
           date[s]['end']=end
           date[s]['sum']=1
           date[s]['e']=e
    xlsxdate(date) #保存到excle


#####################################提取###############################

file_url = '../2-eschColi_B7A-tRNAs.fa.txt'
file_table='date2.xlsx'
date1=[]
date2=[]
date3=[]
i=0


def writedate(date):
    old_excel = xlrd.open_workbook(file_table)  # 打开文件
    #new_excel = copy(old_excel)# 将操作文件对象拷贝，变成可写的workbook对象
    new_workbook = copy(old_excel)
    table =  old_excel.sheet_by_index(0)
    nrows = table.nrows #获取总行数
    ws = new_workbook.get_sheet(0)  # 获取sheet
    for j in date:
        for i in range(1,nrows):
            if table.row(i)[0].value==j[0]:
                ws.write(i, 5, j[12])
                ws.write(i, 6, j[3])
    new_workbook.save(file_table)  # 保存
    pass


with open(file_url, 'r') as f:
    for line in f: #把序列和别的信息放一起
         #print(line)
         i+=1
         for j in line.split():
             date1.append(j)
         if i==3:
             date2.append(date1)
             date1=[]
             i=0
        # print(date)
    for a in range(len(date2)):     #去除文件名字中的< 把序列放一个字符串
        #print(date2[a])
        date2[a][0]=date2[a][0][1:]
        date2[a][len(date2[a])-2]=date2[a][len(date2[a])-2] +date2[a][len(date2[a])-1]
        del date2[a][len(date2[a]) - 1]
        date2[a][3]=date2[a][3][:len(date2[a][3])-1]
        #print(date2[a])
        date3.append(date2[a])
   # print(date3)
    writedate(date3)


#####################################找长度大于等于50###############################

def table1date(date):
    writebook = xlwt.Workbook()  # 打开一个excel
    sheet = writebook.add_sheet('test')  # 在打开的excel中添加一个sheet
    i = 0
    sheet.write(i, 0, '名字')  # 写入名字
    sheet.write(i, 1, '开始位置')  # 开始位置
    sheet.write(i, 2, '结束位置')  # 结束位置
    sheet.write(i, 3, '数量')  # 数量
    sheet.write(i, 4, '正负')
    sheet.write(i, 5, 'trna')  #
    sheet.write(i, 6, 'id')  #
    sheet.write(i, 7, '长度')
    for s in range(len(date)):
        i += 1
        sheet.write(i, 0, date[s][0].value)  # 写入excel，i行0列
        sheet.write(i, 1, date[s][1].value)
        sheet.write(i, 2, date[s][2].value)  # 写入excel，i行0列
        sheet.write(i, 3, date[s][3].value)
        sheet.write(i, 4, date[s][4].value)  # 写入excel，i行0列
        sheet.write(i, 5, date[s][5].value)
        sheet.write(i, 6, date[s][6].value)
        sheet.write(i, 7, int(table.row(i)[2].value) - int(table.row(i)[1].value))
    writebook.save('大肠d2大于等于50.xlsx')  # 一定要记得保存
    print(i)
    pass

def table2date(date):
    writebook = xlwt.Workbook()  # 打开一个excel
    sheet = writebook.add_sheet('test')  # 在打开的excel中添加一个sheet
    i = 0
    sheet.write(i, 0, '名字')  # 写入名字
    sheet.write(i, 1, '开始位置')  # 开始位置
    sheet.write(i, 2, '结束位置')  # 结束位置
    sheet.write(i, 3, '数量')  # 数量
    sheet.write(i, 4, '正负')
    sheet.write(i, 5, 'trna')#
    sheet.write(i, 6, 'id')  #
    sheet.write(i, 7, '长度')
    for s in range(len(date)):
            i += 1
            sheet.write(i, 0, date[s][0].value)  # 写入excel，i行0列
            sheet.write(i, 1, date[s][1].value)
            sheet.write(i, 2, date[s][2].value)  # 写入excel，i行0列
            sheet.write(i, 3, date[s][3].value)
            sheet.write(i, 4, date[s][4].value)  # 写入excel，i行0列
            sheet.write(i, 5, date[s][5].value)
            sheet.write(i, 6, date[s][6].value)
            sheet.write(i, 7, int(table.row(i)[2].value)-int(table.row(i)[1].value))
    writebook.save('大肠d2小于等于50.xlsx')  # 一定要记得保存
    print(i)
    pass

table1=[]
table2=[]

workbook=xlrd.open_workbook('date2.xlsx',on_demand=True)#打开文件
table = workbook.sheet_by_index(0)
nrows = table.nrows  #获取该sheet中的有效行数
for i in range(1,nrows):
    print(table.row(i)[2].value)

    if int(table.row(i)[2].value)-int(table.row(i)[1].value)>=50:  #找数量小于70和大于等于70的
        table1.append(table.row(i))
    else:
        table2.append(table.row(i))

table1date(table1)
table2date(table2)

python计数

该处使用的url网络请求的数据。

---

总结

提示：这里对文章进行总结：
例如：以上就是今天要讲的内容，本文仅仅简单介绍了pandas的使用，而pandas提供了大量能使我们快速便捷地处理数据的函数和方法。