【miRNA-Seq 實戰(zhàn)】二诗祸、準(zhǔn)備分析數(shù)據(jù)跑芳，差異分析以及靶向基因預(yù)測

這里是佳奧！了解了大致流程直颅，我們開始正式的分析實戰(zhàn)博个。

本次實戰(zhàn)的數(shù)據(jù)來源是這篇文章：EBV-miR-BART8-3p induces epithelial-mesenchymal transition and promotes metastasis of nasopharyngeal carcinoma cells through activating NF-κB and Erk1/2 pathways

參考推文：

http://www.bio-info-trainee.com/6997.html

1 miRNA-Seq環(huán)境搭建

參考序列文件

在前一篇我們已經(jīng)準(zhǔn)備好了MIRBASE的miRNA參考序列文件。

mkdir -p ~/refer/miRNA
cd ~/refer/miRNA

wget https://www.mirbase.org/ftp/CURRENT/hairpin.fa.gz ##　38589　reads
wget https://www.mirbase.org/ftp/CURRENT/mature.fa.gz ## 48885 reads 
wget https://www.mirbase.org/ftp/CURRENT/genomes/hsa.gff3

gzip -d hairpin.fa.gz
gzip -d mature.fa.gz

grep sapiens mature.fa |wc -l # 2656 
grep sapiens hairpin.fa |wc # 1917

## Homo sapiens
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
# 這里值得一提的是miRBase數(shù)據(jù)庫下載的序列功偿，居然都是用U表示的盆佣，也就是說就是miRNA序列，而不是轉(zhuǎn)錄成該miRNA的基因序列械荷，而我們測序的都是基因序列共耍。

配置conda

conda create -n mirna
conda activate mirna
conda install -y -c bioconda hisat2 bowtie samtools fastp bowtie2 fastx_toolkit fastqc multiqc

bowtie構(gòu)建索引

##當(dāng)前refer目錄下的文件
(mirna) root 11:02:43 /home/kaoku/refer/miRNA
$ ls -lh
總用量 11M
hairpin.fa
hairpin.human.fa
hsa.gff3
mature.fa
mature.human.fa

bowtie-build mature.human.fa hsa-mature-bowtie-index
bowtie-build hairpin.human.fa hsa-hairpin-bowtie-index

##得到文件
hsa-hairpin-bowtie-index.1.ebwt
hsa-hairpin-bowtie-index.2.ebwt
hsa-hairpin-bowtie-index.3.ebwt
hsa-hairpin-bowtie-index.4.ebwt
hsa-hairpin-bowtie-index.rev.1.ebwt
hsa-hairpin-bowtie-index.rev.2.ebwt
hsa-mature-bowtie-index.1.ebwt
hsa-mature-bowtie-index.2.ebwt
hsa-mature-bowtie-index.3.ebwt
hsa-mature-bowtie-index.4.ebwt
hsa-mature-bowtie-index.rev.1.ebwt
hsa-mature-bowtie-index.rev.2.ebwt

2 文章數(shù)據(jù)下載

下載SRA文件

vdb-config -i 
##進入配置文件后按c,在location of user-repository設(shè)置下載好的文件的位置,設(shè)置完成后按s保存x退出界面

##添加到環(huán)境
export PATH="$PATH:/home/kaoku/biosoft/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin"

##下載第一個項目p1的SRA文件
SRR1542714
SRR1542715
SRR1542716
SRR1542717
SRR1542718
SRR1542719

cat srr.list | while read id; do ( prefetch $id & ); done

mv SRR* /home/kaoku/project/miRNA/p1/

##第二個項目p2的SRA文件
SRR7707744
SRR7707745
SRR7707746
SRR7707747
SRR7707748
SRR7707749
SRR7707750
SRR7707751
SRR7707752
SRR7707753
SRR7707754

cat srr.list | while read id; do ( prefetch $id & ); done

mv SRR* /home/kaoku/project/miRNA/p2/

轉(zhuǎn)fastq文件

##p1目錄
echo {14..19} |sed 's/ /\n/g' |while read id; \
do ( fastq-dump -O ./ --gzip --split-3 SRR15427${id} & ) ; \
done

##p2目錄
echo {44..54} |sed 's/ /\n/g' |while read id; \
do ( fastq-dump -O ./ --gzip --split-3 SRR77077${id} & ) ; \
done

質(zhì)控和清洗測序數(shù)據(jù)

##新建qc目錄，進行質(zhì)量控制
ls ../raw_fq/*gz | xargs fastqc -t 10 -o  ./
multiqc ./

##清洗主要是fastx_toolkit修剪
ls *gz |while read id
do
echo $id 
# fastqc $id 
zcat $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 ${id%%.*}_clean.fq.gz 
done

##會進行兩步操作养葵，質(zhì)量控制并修剪長度
ls *gz |while read id
do
echo $id 
zcat $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 ;
done

#fastq_quality_filter和fastx_trimmer兩個命令征堪，都是來自于fastx_toolkit軟件包：

#fastq_quality_filter 代表根據(jù)質(zhì)量過濾序列
#fastx_trimmer 代表縮短FASTQ或FASTQ文件中的讀數(shù)，根據(jù)質(zhì)量修剪（剪切）序列关拒。

##文件大小變化情況
-rw-r--r-- 1 root root  39M  8月 29 16:06 SRR1542714_clean.fq.gz
-rw-r--r-- 1 root root  50M  8月 29 14:50 SRR1542714.fastq.gz

第一種比對策略：bowtie比對 miRBaes數(shù)據(jù)庫下載的序列 +定量

##同時比對mature與hairpin佃蚜，注意hsa-mature-bowtie-index是文件名前綴，不是目錄名
mature=/home/kaoku/refer/miRNA/hsa-mature-bowtie-index
hairpin=/home/kaoku/refer/miRNA/hsa-hairpin-bowtie-index

ls *_clean.fq.gz |while read id
do
echo $id 
bowtie -n 0 -m1 --best --strata $mature $id -S ${id}_matrue.sam 
bowtie -n 0 -m1 --best --strata $hairpin $id -S ${id}_hairpin.sam 
done

##壓縮.sam成.bam并刪除.sam
ls *.sam|while read id ;do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done
rm *.sam

-n 0 不允許錯配着绊，但是結(jié)果來看比對率很低谐算，修改成允許一個錯配的話：

ls *_clean.fq.gz |while read id
do
echo $id 
bowtie -n 1 -m1 --best --strata $mature $id -S ${id}_matrue.n1.sam 
bowtie -n 1 -m1 --best --strata $hairpin $id -S ${id}_hairpin.n1.sam 
done

隨后得到了.bam文件

(mirna) root 16:46:10 /home/kaoku/project/miRNA/p1/align
$ ls
SRR1542714_clean.fq.gz_hairpin.bam  SRR1542715_clean.fq.gz_matrue.bam   SRR1542717_clean.fq.gz_hairpin.bam  SRR1542718_clean.fq.gz_matrue.bam
SRR1542714_clean.fq.gz_matrue.bam   SRR1542716_clean.fq.gz_hairpin.bam  SRR1542717_clean.fq.gz_matrue.bam   SRR1542719_clean.fq.gz_hairpin.bam
SRR1542715_clean.fq.gz_hairpin.bam  SRR1542716_clean.fq.gz_matrue.bam   SRR1542718_clean.fq.gz_hairpin.bam  SRR1542719_clean.fq.gz_matrue.bam

##查看每個基因表達(dá)量
samtools view .bam | cut -f 3 | sort | uniq -c

批量走 samtools idxstats 得到表達(dá)量矩陣：

ls *.bam | xargs -i samtools index {}
ls *.bam | while read id ;do (samtools idxstats ${id} > ${id}.txt );done
# samtools view matrue.bam |cut -f 3 |sort |uniq -c

第二種比對策略：bowtie2比對參考基因組+定量

下載hg38參考基因組：

wget -c https://genome-idx.s3.amazonaws.com/bt/GRCh38_noalt_as.zip

比對參考基因組

##查看一下index目錄
(mirna) root 19:22:44 /home/kaoku/refer/miRNA/hg38
$ ls
GRCh38_noalt_as.zip  hg38.1.bt2  hg38.2.bt2  hg38.3.bt2  hg38.4.bt2  hg38.rev.1.bt2  hg38.rev.2.bt2

index=/home/kaoku/refer/miRNA/hg38/hg38

ls *_clean.fq.gz |while read id
do
echo $id 
bowtie2 -p 4 -x $index -U $id | samtools sort -@ 4 -o ${id}_genome.bam -
done

##比對結(jié)果
SRR1542714_clean.fq.gz_genome.bam  SRR1542716_clean.fq.gz_genome.bam  SRR1542718_clean.fq.gz_genome.bam
SRR1542715_clean.fq.gz_genome.bam  SRR1542717_clean.fq.gz_genome.bam  SRR1542719_clean.fq.gz_genome.bam

定量分析

##確保安裝subread

gtf=/home/kaoku/refer/miRNA/hsa.gff3

featureCounts -T 4 -F gff -M -t miRNA -g Name -a $gtf -o all.counts.mature.txt *genome* 1>counts.mature.log 2>&1

featureCounts -T 4 -F gff -M -t miRNA_primary_transcript -g Name -a $gtf -o all.counts.hairpin.txt *genome* 1>counts.hairpin.log 2>&1

##結(jié)果
all.counts.hairpin.txt
all.counts.hairpin.txt.summary
all.counts.mature.txt 
all.counts.mature.txt.summary
counts.hairpin.log
counts.mature.log

3 比較兩種定量分析方法差異

比較miRBaes數(shù)據(jù)庫下載的序列和參考基因組的兩組表達(dá)矩陣。

QQ截圖20220829203444.png

hairpin<-read.table('all.counts.hairpin.txt',header = T)
m1<-hairpin[,7:12]
rownames(m1)<-hairpin[,1]

paste *_clean.fq.gz_hairpin.bam.txt | cut -f 1,3,7,11,15,19,23 > hairpin.txt

m2<-read.table('hairpin.txt',header = F,row.names = 1)
ids<-intersect(rownames(m1),rownames(m2))

##構(gòu)建merge matrix
overM<-cbind(m1[ids,],m2[ids,])
##LOG轉(zhuǎn)換后的COUNTS矩陣
pheatmap::pheatmap(cor(log2(overM+1)))
QQ截圖20220829211943.png
##取CPM后的矩陣
cpmM<-log(edgeR::cpm(overM)+1)
pheatmap::pheatmap(cor(cpmM))

##選取TOP500的MAD基因
highM=cpmM[names(tail(sort(apply(cpmM,1,mad)),500)),]
pheatmap::pheatmap(cor(highM))

##獨立選取各自的top500的mad基因归露，再去交集
hg1=names(tail(sort(apply(cpmM[,1:6],1,mad)),500))
hg2=names(tail(sort(apply(cpmM[,7:12],1,mad)),500))
hg=intersect(hg1,hg2)
hgM=cpmM[hg,]
pheatmap::pheatmap(cor(hgM))
QQ截圖20220829212715.png

流程對于定量分析結(jié)果影響很大洲脂。

##進一步探索兩種方法對于比對數(shù)的差異
colnames(overM)=NULL

換讀取另一個

mature<-read.table('all.counts.mature.txt',header = T)
m1<-mature[,7:12]
rownames(m1)<-mature[,1]
cpmM<-log(edgeR::cpm(overM)+1)
pheatmap::pheatmap(cor(cpmM))

m2<-read.table('mature.txt',header = F,row.names = 1)

ids<-intersect(rownames(m1),rownames(m2))
overM<-cbind(m1[ids,],m2[ids,])
pheatmap::pheatmap(cor(log2(overM+1)))

highM=cpmM[names(tail(sort(apply(cpmM,1,mad)),500)),]
pheatmap::pheatmap(cor(highM))

hg1=names(tail(sort(apply(cpmM[,1:6],1,mad)),500))
hg2=names(tail(sort(apply(cpmM[,7:12],1,mad)),500))
hg=intersect(hg1,hg2)
hgM=cpmM[hg,]
pheatmap::pheatmap(cor(hgM))

colnames(overM)=NULL

paste *_clean.fq.gz_matrue.bam.txt | cut -f 1,3,7,11,15,19,23 > mature.txt
QQ截圖20220829213919.png

4 DEG差異分析

##從mature矩陣開始
rm(list = ls())
options(stringsAsFactors = F)
mature<-read.table('all.counts.mature.txt',header = T)
m1<-mature[,7:12]
rownames(m1)<-mature[,1]

##判斷分組
pheatmap::pheatmap(cor(m1))

gl<-rep(c('control','ET1'),3)

##DEG分析
group_list<-gl
exprSet<-m1

suppressMessages(library(DESeq2)) 
(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds)

res <- results(dds, 
               contrast=c("group_list","ET1","control"))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG=as.data.frame(resOrdered)

DEG <- na.omit(DEG)

##判斷上下調(diào)基因
df=DEG
colnames(df)
df$V<- -log10(df$pvalue) 
library(ggpubr)
ggscatter(df, x = "log2FoldChange", y = "V", size=0.5)

df$g<-ifelse(df$pvalue>0.01, 'stable',
             ifelse( df$log2FoldChange > 2 , 'up',
                    ifelse( df$log2FoldChange < - 2, 'down', 'stable')))  
table(df$g)

down stable     up 
 43    492     35

5 靶向基因預(yù)測

參考推文：

https://mp.weixin.qq.com/s/n_UncYeGIQFLneTMK2rTXQ

##載入R包
BiocManager::install("multiMiR",ask = F,update = F)
save(df,file = 'deg.R')

rm(list = ls())
options(stringsAsFactors = F)

load(file = 'deg.R')
library(multiMiR)
db.ver = multimir_dbInfoVersions()
db.ver[,1:3]

下面是使用edgeR包，對普通的轉(zhuǎn)錄組counts表達(dá)矩陣（miRNA）做差異分析剧包，并且拿到感興趣的miRNA基因集：

library(edgeR)
library(multiMiR)

# Load data
counts_file  <- system.file("extdata", "counts_table.Rds", package = "multiMiR")
strains_file <- system.file("extdata", "strains_factor.Rds", package = "multiMiR")
counts_table   <- readRDS(counts_file)
strains_factor <- readRDS(strains_file)
table(strains_factor)

# Standard edgeR differential expression analysis
design <- model.matrix(~ strains_factor)

# Using trended dispersions
dge <- DGEList(counts = counts_table)
dge <- calcNormFactors(dge)
dge$samples$strains <- strains_factor
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)

# Fit GLM model for strain effect
fit <- glmFit(dge, design)
lrt <- glmLRT(fit)

# Table of unadjusted p-values (PValue) and FDR values
p_val_DE_edgeR <- topTags(lrt, adjust.method = 'BH', n = Inf)

# Getting top differentially expressed miRNA's
top_miRNAs <- rownames(p_val_DE_edgeR$table)[1:10]

有了感興趣的miRNA基因集恐锦，就可以查詢它們的靶基因

library(multiMiR)
# Plug miRNA's into multiMiR and getting validated targets
multimir_results <- get_multimir(org     = 'mmu',
                                 mirna   = top_miRNAs,
                                 table   = 'validated',
                                 summary = TRUE)
head(multimir_results@data)
table(multimir_results@data$mature_mirna_id)
dim(multimir_results@data)

相關(guān)的內(nèi)容就到這里，miRNA-Seq 實戰(zhàn)課程就到這里疆液。

我們下一個篇章（單細(xì)胞篇）再見Ｒ磺Α！堕油！

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末潘飘，一起剝皮案震驚了整個濱河市肮之，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌卜录，老刑警劉巖戈擒，帶你破解...
    沈念sama閱讀 218,525評論 6贊 507
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異艰毒，居然都是意外死亡筐高，警方通過查閱死者的電腦和手機，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,203評論 3贊 395
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門现喳，熙熙樓的掌柜王于貴愁眉苦臉地迎上來凯傲，“玉大人犬辰，你說我怎么就攤上這事嗦篱。” “怎么了幌缝？”我有些...
    開封第一講書人閱讀 164,862評論 0贊 354
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵灸促，是天一觀的道長。 經(jīng)常有香客問我涵卵，道長浴栽，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 58,728評論 1贊 294
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任轿偎，我火速辦了婚禮典鸡，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘坏晦。我一直安慰自己萝玷，他們只是感情好，可當(dāng)我...
    茶點故事閱讀 67,743評論 6贊 392
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布昆婿。 她就那樣靜靜地躺著球碉，像睡著了一般。 火紅的嫁衣襯著肌膚如雪仓蛆。 梳的紋絲不亂的頭發(fā)上睁冬，一...
    開封第一講書人閱讀 51,590評論 1贊 305
  • 城市分裂傳說
    那天，我揣著相機與錄音看疙，去河邊找鬼豆拨。 笑死，一個胖子當(dāng)著我的面吹牛能庆，可吹牛的內(nèi)容都是我干的施禾。 我是一名探鬼主播，決...
    沈念sama閱讀 40,330評論 3贊 418
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼相味，長吁一口氣：“原來是場噩夢啊……” “哼拾积！你這毒婦竟也來了殉挽？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 39,244評論 0贊 276
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤拓巧，失蹤者是張志新（化名）和其女友劉穎斯碌，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體肛度，經(jīng)...
    沈念sama閱讀 45,693評論 1贊 314
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡傻唾，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,885評論 3贊 336
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了承耿。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片冠骄。...
    茶點故事閱讀 40,001評論 1贊 348
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖加袋，靈堂內(nèi)的尸體忽然破棺而出凛辣，到底是詐尸還是另有隱情，我是刑警寧澤职烧，帶...
    沈念sama閱讀 35,723評論 5贊 346
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布扁誓，位于F島的核電站，受9級特大地震影響蚀之，放射性物質(zhì)發(fā)生泄漏蝗敢。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點故事閱讀 41,343評論 3贊 330
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一足删、第九天 我趴在偏房一處隱蔽的房頂上張望寿谴。 院中可真熱鬧，春花似錦失受、人聲如沸讶泰。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,919評論 0贊 22
  • 一樁弒父案贱纠，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽峻厚。三九已至，卻和暖如春谆焊，著一層夾襖步出監(jiān)牢的瞬間惠桃，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,042評論 1贊 270
  • 情欲美人皮
    我被黑心中介騙來泰國打工辖试， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留辜王，地道東北人。 一個月前我還...
    沈念sama閱讀 48,191評論 3贊 370
  • 代替公主和親
    正文 我出身青樓罐孝，卻偏偏與公主長得像呐馆，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子莲兢，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 44,955評論 2贊 355

推薦閱讀更多精彩內(nèi)容