一般流程：参考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%左右。
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档
前言
? ? ? ?随着测序技术的发展，数据分析这门技术也越来越重要，作为生信小白自学着实有些艰难，分享一篇学习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 <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>
实例
$ 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
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 <annotation_file> -o <output_file> 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提供了大量能使我们快速便捷地处理数据的函数和方法。