深夜免费网站 miRNA测序数据的上游定量历程实战演练

最近miRNA范畴果然取得了诺贝尔奖深夜免费网站，是以我也凑一下骚扰念念学习，是以就委用神通高大的曾赤诚提供了一个miRNA测序数据案例，这个数据集是GSE181922，其中一共包括了40例样品的的miRNA数据。如下所示：

图片

miRNA的上游分析历程跟mRNA的上游历程很相似：环境部署——数据下载——稽察数据(非质控)——数据质控清洗——数据比对——数据定量https://www.bilibili.com/video/BV1zK411n7qr 

图片

1.基于conda的环境部署/软件装配：

尝试使用ARM架构(M1/M2芯片) 去装配fastqc trim-galore hisat2 subread multiqc samtools salmon fastp，发现这些软件中有几个是不兼容的。是以需要改回底本的x86_64架构(Intel芯片)，如果非mac/M1/M2的不需要用这种神态。

CONDA_SUBDIR=osx-64 conda create -n miRNA_x86_64 python=3.9  conda activate miRNA_x86_64  conda install -y -c bioconda sra-tools hisat2 bowtie samtools fastp bowtie2 fastx_toolkit fastqc  conda install -y -c hcc aspera-cli  conda install -y sra-tools  

2.下载相应数据库数据

miRbase是miRNA推敲范畴内最泰斗的数据库之一，提供了miRNAs序列以及防卫，定位，发夹序列等信息，以及提供定名做事。

mkdir mirna  mkdir reference  cd ./reference  #nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa &  #nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa &  #gunzip hairpin.fa.gz  #gunzip mature.fa.gz    #不知为何，笔者这边一直出现聚首失败，委果没见识就奏凯从官网进行了点击下载

1.  前体 miRNA（hairpin.fa）：

识别新 miRNA： 通过比对发夹状序列，推敲东谈主员不错瞻望或识别新的 miRNA，因为腾达成的 miRNA 在细胞内领先造成发夹结构。结构分析： 发夹状结构是 miRNA 独有的二级结构，通过分析它的结构和序列特征，不错更好地了解 miRNA 的生成机制。料理和调动为锻练 miRNA： 前体 miRNA 是锻练 miRNA 的来源，前体文献不错匡助模拟或推敲 miRNA 在细胞内的生成过程，举例 Dicer 酶切割的具体位置和生成机制。

2.  锻练 miRNA（mature.fa）：

功能性分析： 锻练 miRNA 是奏凯调控基因抒发的功能分子，不错招引到特定的 mRNA 靶位点。mature.fa 文献不错用于 miRNA 靶向基因瞻望、通路分析和功能推敲。抒发定量： 在内容的抒发定量分析（如 RNA-seq）中，比对锻练 miRNA 序列不错匡助准确识别 miRNA，并进行定量，从而用于下流的各异抒发分析。基因调控网罗： 锻练 miRNA 文献可用于构建 miRNA-mRNA 调控网罗，推敲 miRNA 在特定生物学过程中的作用。

在这里这两个文献的作用主如果进行序列比对。

图片

3.Check 下载到腹地的数据

通达hairpin.fa文献不错看到数据的时期

cel-let-7 是序列称号。MI0000001 是 miRBase 数据库中对应的独一 ID。Caenorhabditis elegans 是该序列的来源物种。let-7 stem-loop 形色了该序列是 let-7 miRNA 的发夹环结构。紧接的两行是 let-7 的核苷酸序列。

图片

cat hairpin.fa | grep '>'| awk '{print $3,$4}'| sort |uniq -c | wc -l  # 265

cat hairpin.fa: 读取 hairpin.fa 文献的一齐内容，并输出到终局。grep '>': 筛选出包含 > 字符的行。FASTA 时期中，> 起原的行默示序列的防卫信息，如 miRNA 称号和其他信息，而不是序列自身。awk '{print 4}': 这里 {print

图片

4} 默示输出第三个和第四个字段（即以空格或制表符分隔的第三和第四部分）。在 miRNA FASTA 文献中，第三个和第四个字段可能是与 miRNA 称号和种类干系的信息。sort对提真金不怕火的第三和第四字段进行排序。uniq -c: 统计每个独一的第三、四字段组合的出现次数。uniq -c 会对换取的行进行计数。举例，如果 miRNA_type1 出现了屡次，则会输出雷同 5 miRNA_type1。在使用 uniq 之前，必须先对内容进行 sort，不然无法识别换取的行。wc -l：统计输出的行数。wc -l 统计 uniq -c 输出的总行数，即不同 miRNA 类型组合的数目。

cat mature.fa | grep '>'| awk '{print $3,$4}'| sort |uniq -c | wc -l  # 265

接着不雅察东谈主类这个物种的miRNA的数目

grep sapiens mature.fa |wc -l  # 2656  

grep sapiens mature.fa | wc -l：grep sapiens mature.fa：从文献 mature.fa 中查找包含 "sapiens" 的行。| wc -l：将匹配的行数通过管谈传递给 wc -l，统计这些行的数目。最终复返包含 "sapiens" 的总行数。grep sapiens hairpin.fa | wc：grep sapiens hairpin.fa：从文献 hairpin.fa 中查找包含 "sapiens" 的行。| wc：wc 会复返三个值：行数、单词数、字节数。由于莫得加 -l 参数，成果中会包含所有这三个统计值，列出包含 "sapiens" 的行数、单词总和以及字符总和。

接着不雅察有若干序列，4行为一条序列

zless -S mature.fa | paste - - - - |wc -l  # 24443  zless -S hairpin.fa | paste - - - - |wc -l  # 30029

接着查验一下前体和锻练体长度：

前体miRNA和锻练体miRNA：前体miRNA长度一般是70-120碱基，时常是茎环(发夹，hairpin)结构。锻练之后一般是22个碱基。(曾赤诚的perl的示例代码)

# 前体长度  awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < hairpin.fa | tr "\t" "\n" | grep -v '>' | awk '{print length}' | uniq -c | sort -n -k2    # 锻练体长度  awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < mature.fa | tr "\t" "\n" | grep -v '>' | awk '{print length}' | uniq -c | sort -n -k2

4.构建索引

构建 miRNA 序列的索引库（举例使用 bowtie 构建 hairpin.fa 和 mature.fa 的索引）不错显赫升迁后续比对过程的速率和准确性，比如：1. 加速比对过程；2. 减少诡计资源需求；3. 升迁比瞄准确性；

U->T退换

为什么要进行U-> T退换：在 RNA 序列中，碱基用“U”（尿嘧啶）默示，而 DNA 序列中对应的是“T”（胸腺嘧啶）。大大都比对器具，如 Bowtie，主如果针对 DNA 序列联想的，因此它们默意识别“ATCG”四种碱基。在这种情况下，如果不将 RNA 中的“U”退换为“T”，比对器具会无法正确识别和比对 RNA 序列。

perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print}' hairpin.fa > hairpin.human.fa  perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print}' mature.fa > mature.human.fa

perl -alne '...'： perl：调用 Perl 剧本。-a：启用自动分段形状，将每行分割成字段，保存在 @F 数组中（这里未用到 @F）。-l：自动料理换行符，使输出更整王人。-n：轮回读取每一转，但不会自动打印输出。if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};：/^>/ ：检测行是否以 > 起原，这时常默示 FASTA 时期中的序列 ID 行。**if(/Homo/)**：如果 ID 行中包含“Homo”（指代东谈主类干系的序列），则将 $tmp 确立为 1（即允许输出该序列的标记）；不然确立为 0（放弃该序列）。这一步细目每条序列是否为东谈主类序列，仅料理包含 Homo 的序列。next if $tmp!=1;： next if ：如果tmp 不是 1，跳过该行，即非东谈主类序列奏凯跳过不意理。s/U/T/g if !/>/;：s/U/T/g：将行中的“U”替换为“T”，全局替换（即一转中所有“U”均替换为“T”）。if !/>/：仅当行中不含 > 绚烂时践诺替换，即欺诈在内容序列行，而非序列 ID 行。print：print：打印料理后的行。对相宜条件的序列和序列 ID 均输出到指定文献。> hairpin.human.fa 和 > mature.human.fa：>将圭臬输出重定向到文献，辞别保存到 hairpin.human.fa 和 mature.human.fa。

图片

Bowtie和Bowtie2的区别是什么：

Bowtie：领受基于 Burrows-Wheeler 变换（BWT）和 FM-index 的算法，适当对短序列（时常为小于 50bp 的短 RNA 或短读长 DNA）进行快速比对。Bowtie2：领受了更复杂的比对算法，使用 Burrows-Wheeler 变换和动态诡沟通法来救助长片断的局部和全局比对，因此适当较长的读长（一般在 50bp 以上），包括 DNA 和 RNA-seq 数据。

bowtie-build hairpin.human.fa hsa-hairpin-bowtie-index  bowtie-build mature.human.fa hsa-mature-bowtie-index    # check  ls -lh

图片

5.下载数据

勾选念念要下载的数据，并点击accession list，并把list放入mirna文献夹中

图片

cd ../  cd ./mirna    # check  ls -lh

使用prefetech下载数据，这里需要把SRRlist和SRA toolkit软件装配好。除了这种神态，咱们也不错遴荐aspera下载神态

nohup bash -c 'cat SRR_Acc_List.txt | while read id; do  echo "Downloading $id"  prefetch $id  done' &> download.log &  

把sra数据批量退换为fastq数据

# 领先需要知谈fastq-dump器具的位置  which fastq-dump  # /opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump    # 稽察文献夹中的数据是如何样的  ls    # 要明确一下echo和SRR的ID情况  # 输入进终局的时分一定要再三明确代码情况  dump=/opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump  echo {02..25} |sed 's/ /\n/g' |while read id; \  do ( $dump -O ./ --gzip --split-3 SRR154179${id}) ;\  done    # 数据有点大 笔者就分开下载了  dump=/opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump  echo {35..50} | sed 's/ /\n/g' | while read id; do  ($dump -O ./ --gzip --split-3 SRR154179${id})  done

{02..25} 会生成一个从 02 到 25 的数字序列。sed 's/ /\n/g' 将生成的序列号中的空格替换为换行符，以便逐行读取数字。while read id; do ... done 造成一个轮回，逐行读取序列号并存储在变量 id 中。其中 {id}：指定待退换的 SRA 文献，${id} 会插入轮回读取的数字部分，生成雷同 SRR15417902、SRR15417903 等文献名。

图片

6.数据质控和清洗

数据质控稽察

# 对面前文献夹中所有以fastq.gz文献进行质料规定  fastqc -t 2 -o ./ *.fastq.gz  # 对面前文献夹中所有以fastq.gz文献进行整合质料规定，输出到fastq_qc文献夹中  multiqc ./*zip -o ./fastq_qc

fastqc：调用 FastQC 器具，用于分析测序数据的质料。-t 2：指定使用 2 个线程并行料理文献，以加速分析速率。-o ./：指定输出目次为面前目次（./），FastQC 生成的发挥文献将保存在面前目次中。*.fastq.gz：匹配面前目次下所有以 .fastq.gz 落幕的文献，手脚输入文献进行质料规定分析。

图片

认真数据清洗

# 查验文献压缩时期  file /Users/zaneflying/Desktop/miRNA_Z/mirna/SRR15417902.fastq.gz    # trim+clean  # 著作用了miRquant 2.0这个器具进行trim，但笔者已经按照曾赤诚的历程进行料理  ls /Users/zaneflying/Desktop/miRNA_Z/mirna/*.gz | while read id; do  echo $id  gzip -cd $id | fastq_quality_filter -v -q 20 -p 80 -Q33 -i - -o tmp  fastx_trimmer -v -f 1 -l 27 -m 15 -i tmp -Q33 -z -o ${id%%.*}_clean.fq.gz  fastqc -t 2 -o ./ ${id%%.*}_clean.fq.gz  done  # check一下  ls -lh *_clean.fq.gz

图片

7.数据比对

把柄miRBase数据库下载的序列进行比对和定量。

# mature清洗和定量  index=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa-mature-bowtie-index    ls *_clean.fq.gz | while read id; do  echo $id  bowtie -p 2 -x $index $id -S tmp  samtools view -bS -@ 2 tmp -o ${id}_mature.bam  done    # hairpin清洗和定量  index=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa-hairpin-bowtie-index.    ls *_clean.fq.gz | while read id; do  echo $id  bowtie -p 2 -x $index $id -S tmp  samtools view -bS -@ 2 tmp -o ${id}_hairpin.bam  done

*ls _clean.fq.gz: 列出所有以 _clean.fq.gz 落幕的文献，即预感理过的 miRNA 序列文献。while read id; do ... done: 使用 while 轮回逐一读取文献名并将其赋值给变量 id，然后对每个文献践诺轮回内的敕令。echo $id: 打印面前正在料理的文献名，用于查验程度。bowtie -p 2 -x id -S tmp: -p 2：指定使用 2 个线程来加速料理。-x ：指定的索引文献，index 应该是您之前创建的 miRNA 参考序列的索引（在您的例子中应该是 /Users/zaneflying/Desktop/miRNA_Z/reference/mature）。$id：输入文献，即面前的 _clean.fq.gz 文献。-S tmp**：指定输出文献 tmp，这将是一个 SAM 时期的中间文献。samtools view -bS -@ 2 tmp -o {id}_mature.bam：输出文献，将 BAM 文献定名为输入文献名加 _mature.bam 后缀。

图片

对比成果中发现只消1507条reads对应上，也便是说确实所有都莫得比对上的情况，很蒙胧。应该是我莫得学好：

然后尝试更换一下参考基因组，著作中提到的是hg19

图片

笔者这里使用GRCh38进行对比，不外这个并不是要点哈。下载流GRCh38程可见转录组上游分析历程推文。

图片

# 生成索引  bowtie2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38.dna  # index索引要求  index=/Users/zaneflying/Desktop/miRNA_Z/GRCh38.113/GRCh38.dna  # bowtie运行调动  ls *_clean.fq.gz | while read id  do  echo $id  bowtie -p 2 -x $index $id -S tmp ;  samtools view -bS -@ 2 tmp -o ${id}_genome.bam  done

ls *_clean.fq.gz | while read id：列出面前目次下所有文献名以 _clean.fq.gz 落幕的文献。while read id 用来逐行读取这些文献名，并将文献名存储在变量 id 中。echo $id：打印面前正在料理的文献名，以便跟踪程度。bowtie -p 2 -x id -S tmp：使用 bowtie 对文献 进行基因组比对。指定使用个线程。index 指定索引文献旅途，即 。id 是输入的 FASTQ 文献。-S tmp 将比对成果输出为 SAM 时期，临时保存在文献 tmp 中。samtools view -bS -@ 2 tmp -o {id}_genome.bam 指定输出文献名为 ${id}_genome.bam，即与输入文献同名，但以 _genome.bam 落幕。

图片

这个对比成果情况就对付能“让东谈主接管”啦~

8.数据定量

著作顶用的是miRquant 2.0

图片

笔者使用featurecounts去定量， 需要先去miRBase潦倒载hsa.gff3

图片

# 下载hsa.gff文献,放到reference文献夹中  nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3 &  # 如果网罗不好，就奏凯办动下载    # 装配subread  conda install -c bioconda subread    # 运行比对  gtf=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa.gff3  featureCounts -T 2 -F gff -M -t miRNA -g Name -a $gtf -o all.counts.mature.txt *genome* 1>counts.mature.log 2>&1  featureCounts -T 2 -F gff -M -t miRNA_primary_transcript -g Name -a $gtf -o all.counts.hairpin.txt *genome* 1>counts.hairpin.log 2>&1  

featureCounts：调用 featureCounts 器具进行基因计数。-T 2：指定使用 2 个线程。-F gff：指定输入文献的防卫时期为 GFF。-M：允很多重比对的 reads（即与多个位置比对的 reads）。-t miRNA：在 GFF 文献中，遴荐类型为 miRNA 的要求，这样不错仅对锻练 miRNA 计数。-g Name：遴荐 GFF 文献中 Name 字段手脚基因 ID。-a ：指定防卫文献旅途，即前边确立的gtf 变量。-o all.counts.mature.txt：输出文献称号，包含锻练 miRNA 的计数。genome：指定输入的 BAM 文献，genome 是通配符，默示所知称号包含“genome”的 BAM 文献。1>counts.mature.log 2>&1：将圭臬输出和圭臬诞妄重定向到日记文献 counts.mature.log。第二行敕令与第一转基本换取，独一的区别是 -t 选项为 miRNA_primary_transcript，用于遴荐前体 miRNA 要求，从而统计前体 miRNA 的计数。输出文献为 all.counts.hairpin.txt，并将日记信息输出到 counts.hairpin.log。

图片

因为比对有问题，定量也很难保证，是以拿到了矩阵也很难进行后续分析：

图片

后续的分析基本上等同于转录组测序抒发量矩阵，便是各异分析等统计可视化：

图片

参考良友：Multiomic analysis of microRNA-mediated regulation reveals a proliferative axis involving miR-10b in fibrolamellar carcinoma. JCI Insight. 2022 Jun 8;7(11):e154743.DNAJB1-PRKACA fusion protein-regulated LINC00473 promotes tumor growth and alters mitochondrial fitness in fibrolamellar carcinoma. PLoS Genet 2024 Mar;20(3):e1011216.Chemical, Molecular, and Single-nucleus Analysis Reveal Chondroitin Sulfate Proteoglycan Aberrancy in Fibrolamellar Carcinoma. Cancer Res Commun 2022 Jul;2(7):663-678.生信妙技树：跋文

我如实是看完毕素养视频，以及配套的条记，然而不知谈为什么成果就大相径庭，一个东谈主学习生信便是如斯的无聊和秘密！

小RNA建库测序后的数据分析-实例老师随着生信妙技树学习microRNA测序学习构建miRNA-seq数据分析环境miRNAseq数据分析这样多年了它的历程也莫得固定这5个miRNA构成的肺鳞癌会诊基因集在tcga数据库能复现吗什么，给你了你这样多miRNA靶基因查询R包和网页器具你果然不知谈如何使用对miRNA进行go和kegg等功能数据库数据库防卫使用miRNAtap数据源提真金不怕火miRNA的瞻望靶基因成果谁说Windows下无法作念生信分析（植物miRNA gene瞻望给你看）你但愿这个探针防卫到卵白编码基因已经miRNA的基因呢如果miRNA的3p和5p功能不同样miRNA、LncRNA、CircRNA靠谱小结诡计MiRNA–mRNA抒发干系性使用多个网页器具瞻望MiRNA–mRNA相互作用一篇著作学会miRNA-seq分析

致谢：感谢曾赤诚以及生信妙技树团队整体成员。

注：若对内容有猜忌八成有发现明确诞妄的一又友深夜免费网站，请干系后台(迎接交流)。 

本站仅提供存储做事，所有内容均由用户发布，如发现存害或侵权内容，请点击举报。