<source id="wfbuo"><optgroup id="wfbuo"></optgroup></source>
<rt id="wfbuo"></rt><strong id="wfbuo"></strong>
<rt id="wfbuo"><optgroup id="wfbuo"></optgroup></rt>
<rt id="wfbuo"></rt>

  • <strong id="wfbuo"></strong>

    <menuitem id="wfbuo"></menuitem>
  • <strong id="wfbuo"></strong>
  • 當前位置 博文首頁 > TOP生物信息:一文學會常規轉錄組分析

      TOP生物信息:一文學會常規轉錄組分析

      作者:TOP生物信息 時間:2021-02-20 16:35

      數據來源:Cytosolic acetyl-CoA promotes histone acetylation predominantly at H3K27 in Arabidopsis;GSE79524

      我只試做了轉錄組分析那一部分。簡單概括就是為了評估乙;瘜虮磉_的影響,對野生型和突變體都做了轉錄組分析(基因差異表達分析和GO注釋)。

      原文方法,我換成了自己較熟悉的幾個工具

      1. 數據獲取及質控

      #1.腳本查看
      $ cat dir_6.txt 
      SRR3286802
      SRR3286803
      SRR3286804
      SRR3286805
      SRR3286806
      SRR3286807
      
      $ cat 1.sh 
      dir=$(cut -f1 dir_6.txt)
      for sample in $dir
      do
      prefetch $sample
      done
      
      $ cat 2.sh 
      for i in `seq 2 7`
      do
      fastq-dump --gzip --split-3 --defline-qual '+' --defline-seq '@$ac-$si/$ri' ~/ncbi/public/sra/SRR328680${i}.sra -O /ifs1/Grp3/huangsiyuan/learn_rnaseq/mRNA-seq/data/ &
      done
      
      #2.下載及解壓
      sh 1.sh
      sh 2.sh
      
      #解壓之后是這樣的,可以看出是雙端測序
      $ ls
      1.sh       SRR3286802_1.fastq.gz  SRR3286803_2.fastq.gz  SRR3286805_1.fastq.gz  SRR3286806_2.fastq.gz
      2.sh       SRR3286802_2.fastq.gz  SRR3286804_1.fastq.gz  SRR3286805_2.fastq.gz  SRR3286807_1.fastq.gz
      dir_6.txt  SRR3286803_1.fastq.gz  SRR3286804_2.fastq.gz  SRR3286806_1.fastq.gz  SRR3286807_2.fastq.gz
      
      #3.質量檢測:fastqc,multiqc
      ls *fastq.gz | while read id
      do
      fastqc ${id} &
      done
      multiqc *fastqc*
      
      #4.過濾接頭序列及低質量reads片段
      ##就這組數據來說,質量檢測的結果表明數據質量很好,因此這里省略的過濾這一步。
      

      2. 下載gff/gtf注釋文件并提取出感興趣的基因/轉錄本區間

      一個基因可能對應不同的轉錄本,不同的轉錄本可能對應不同的功能。
      以擬南芥的gff注釋文件為例:

      #提取基因
      $ less Arabidopsis_thaliana.TAIR10.42.gff3 | awk '{ if($3=="gene") print $0 }' | cut -d ";" -f1 | head -n 5
      1	araport11	gene	3631	5899	.	+	.	ID=gene:AT1G01010
      1	araport11	gene	6788	9130	.	-	.	ID=gene:AT1G01020
      1	araport11	gene	11649	13714	.	-	.	ID=gene:AT1G01030
      1	araport11	gene	23121	31227	.	+	.	ID=gene:AT1G01040
      1	araport11	gene	31170	33171	.	-	.	ID=gene:AT1G01050
      
      #提取轉錄本
      $ less Arabidopsis_thaliana.TAIR10.42.gff3 | awk '{ if($3=="mRNA") print $0 }' | cut -d ";" -f1 | head -n 10
      1	araport11	mRNA	3631	5899	.	+	.	ID=transcript:AT1G01010.1
      1	araport11	mRNA	6788	8737	.	-	.	ID=transcript:AT1G01020.6
      1	araport11	mRNA	6788	8737	.	-	.	ID=transcript:AT1G01020.2
      1	araport11	mRNA	6788	9130	.	-	.	ID=transcript:AT1G01020.3
      1	araport11	mRNA	6788	9130	.	-	.	ID=transcript:AT1G01020.5
      1	araport11	mRNA	6788	9130	.	-	.	ID=transcript:AT1G01020.4
      1	araport11	mRNA	6788	9130	.	-	.	ID=transcript:AT1G01020.1
      1	araport11	mRNA	11649	13714	.	-	.	ID=transcript:AT1G01030.1
      1	araport11	mRNA	11649	13714	.	-	.	ID=transcript:AT1G01030.2
      1	araport11	mRNA	23121	31227	.	+	.	ID=transcript:AT1G01040.1
      

      從ID可以看出AT1G01020基因有6個轉錄本。這篇文獻比較的是基因層面的表達差異,所以這里我提取的是基因gff,算上線粒體和葉綠體上的基因一共27655個。

      $ less Arabidopsis_thaliana.TAIR10.42.gff3 | awk '{ if($3=="gene") print $0 }' > gene27655.gff
      

      3. 將reads比對到參考基因組

      3.1 為參考基因組建立索引
      nohup ~/mysoft/hisat2-2.1.0/hisat2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Arabidopsis_thaliana &
      

      3.2 比對
      $ cat 3.sh 
      for i in `seq 2 7`
      do
      hisat2 -x ~/learn_rnaseq/mRNA-seq/ref/Arabidopsis_thaliana -p 8 \
      -1 ~/learn_rnaseq/mRNA-seq/data/SRR328680${i}_1.fastq.gz \
      -2 ~/learn_rnaseq/mRNA-seq/data/SRR328680${i}_2.fastq.gz \
      -S ~/learn_rnaseq/mRNA-seq/map/SRR328680${i}.sam
      done
      
      $ nohup sh 3.sh &
      

      從報告文件來看比對率都挺高的,97%以上。

      3.3 sam轉bam, 并排序
      $ cat 1.sh 
      for i in `seq 2 7`
      do
      /ifs1/Software/biosoft/samtools-1.9/samtools view -@ 8 -Sb SRR328680${i}.sam > SRR328680${i}.bam
      /ifs1/Software/biosoft/samtools-1.9/samtools sort -@ 8 -n SRR328680${i}.bam > SRR328680${i}.n.bam
      done
      $ nohup  sh 1.sh &
      

      4. 計算表達量

      Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ...
      nohup ~/mysoft/subread-1.6.3-source/bin/featureCounts -F GTF -t gene -g gene_id -T 4 -a ~/learn_rnaseq/mRNA-seq/ref_ht/gene27655.gff -o ~/learn_rnaseq/mRNA-seq/count/test ~/learn_rnaseq/mRNA-seq/map/SRR328680*.n.bam &
      #上面這一步運行完之后會多出test.summary,test這兩個文件
      nohup ~/mysoft/subread-1.6.3-source/bin/featureCounts -F GTF -t gene -g gene_id -T 4 -p -a ~/learn_rnaseq/mRNA-seq/ref_ht/gene27655.gff -o ~/learn_rnaseq/mRNA-seq/count/test2 ~/learn_rnaseq/mRNA-seq/map/SRR328680*.n.bam 2> log2 &
      #上面這行只多了-p選項,為了看看在雙端測序中統計fragments和reads有什么區別。運行完了多出test2.summary,test2兩個文件。
      

      從test2和text兩個矩陣文件中,可以看出加了p之后求的是fragments的計數,數值大概是reads計數的1/2,這很好理解,因為fragment的定義中就是包含了一對reads的。

      $ lsx test2 |head -n 20 | tail -n 5
      AT1G01140       1       64166   67774   -       3609    2309    2073    2323    1522    1742    1254
      AT1G01150       1       69911   72138   -       2228    3       0       0       1       2       4
      AT1G01160       1       72339   74096   +       1758    766     730     857     745     819     747
      AT1G01170       1       73931   74737   -       807     1796    1758    1748    1061    1540    1340
      AT1G01180       1       75390   76845   +       1456    374     315     311     358     284     312
      
      $ lsx test |head -n 20 | tail -n 5
      AT1G01140       1       64166   67774   -       3609    4587    4118    4613    3025    3469    2484
      AT1G01150       1       69911   72138   -       2228    6       0       0       2       4       8
      AT1G01160       1       72339   74096   +       1758    1449    1383    1637    1412    1538    1409
      AT1G01170       1       73931   74737   -       807     3220    3121    3169    1914    2771    2369
      AT1G01180       1       75390   76845   +       1456    745     619     619     708     558     621
      

      將test2文件中的這7列提取出來即為表達矩陣。

      grep "^AT" test2 | cut -f1,7,8,9,10,11,12 > 6sample_count_matrix.txt
      

      5. 差異表達分析

      6sample_count_matrix.txt

      使用DESeq2 包。

      a <- as.matrix(read.table("6sample_count_matrix.txt",sep="\t",header = T,row.names=1))
      condition <- factor(c(rep("WT",3),rep("acc1.5",3)), levels = c("WT","acc1.5"))
      colData <- data.frame(row.names=colnames(a), condition)
      #查看一下
      > colData
              condition
      WT1            WT
      WT2            WT
      WT3            WT
      acc1.51    acc1.5
      acc1.52    acc1.5
      acc1.53    acc1.5
      
      #確保
      > all(rownames(colData) == colnames(a))
      [1] TRUE
      
      dds <- DESeqDataSetFromMatrix(a, colData, design= ~ condition)
      dds <- DESeq(dds)
      #運行結束會報告
      estimating size factors
      estimating dispersions
      gene-wise dispersion estimates
      mean-dispersion relationship
      final dispersion estimates
      fitting model and testing
      
      #得到初步結果并查看
      res <-  results(dds, contrast=c("condition", "acc1.5", "WT"))
      ##log2(fold change)也是一個衡量差異顯著性的指標
      resLFC <- lfcShrink(dds, coef="condition_acc1.5_vs_WT", type="apeglm")
      resLFC
      ##按照p值由小到大排列
      resOrdered <- res[order(res$pvalue),]
      resOrdered
      ##矯正p值小于0.001的有多少個差異基因
      sum(res$padj < 0.001, na.rm=TRUE)
      ##畫個MA-plot看看大體趨勢
      plotMA(res, ylim=c(-2,2))
      plotMA(resLFC, ylim=c(-2,2))
      
      #按照自定義的閾值提取差異基因并導出
      diff_gene_deseq2 <-subset(res, padj < 0.001 & abs(log2FoldChange) > 1)
      write.csv(diff_gene_deseq2,file= "DEG_acc1.5_vs_WT.csv")
      

      6. 簡單的GO注釋

      首選clusterProfiler包。

      library(”clusterProfiler“)
      library("org.At.tair.db")
      e <- read.table("DEG_acc1.5_vs_WT.csv", header = T,sep = ",")
      f <- e[,1]
      #轉換ID并將ENTREZID編號作為后面的識別ID
      ids <- bitr(f, fromType="TAIR", toType=c("SYMBOL","ENTREZID"), OrgDb="org.At.tair.db")
      f <- ids[,3]
      
      按照GO系統給基因分類
      ggo <- groupGO(gene     = f,
                     OrgDb    = org.At.tair.db,
                     ont      = "CC",
                     level    = 3,
                     readable = TRUE)
      > head(ggo)
                         ID                    Description Count GeneRatio
      GO:0005886 GO:0005886                plasma membrane   237   237/718
      GO:0005628 GO:0005628              prospore membrane     0     0/718
      GO:0005789 GO:0005789 endoplasmic reticulum membrane     3     3/718
      GO:0019867 GO:0019867                 outer membrane     6     6/718
      GO:0031090 GO:0031090             organelle membrane    63    63/718
      GO:0034357 GO:0034357        photosynthetic membrane    23    23/718
      
      富集分析
      ego2 <- enrichGO(gene         = ids$TAIR,
                       OrgDb         = org.At.tair.db,
                       keyType       = 'TAIR',
                       ont           = "CC",
                       pAdjustMethod = "BH",
                       pvalueCutoff  = 0.01,
                       qvalueCutoff  = 0.05,
                       readable      = TRUE)
      
      > head(ego2)
                         ID          Description GeneRatio   BgRatio       pvalue    p.adjust      qvalue
      GO:0009505 GO:0009505 plant-type cell wall    20/654 274/24998 3.989815e-05 0.002907596 0.002713546
      GO:0048046 GO:0048046             apoplast    22/654 328/24998 5.995043e-05 0.002907596 0.002713546
                                                                                                                                                   geneID
      GO:0009505         ATBXL2/ATPMEPCRA/LRX1/AtWAK1/ATPME2/GLL22/LRX2/SS3/ATPAP1/ATC4H/CASP1/BGLU15/RHS12/BGAL1/ATGRP-5/ATPCB/EARLI1/ATPGIP2/FLA13/PME5
      GO:0048046 iPGAM1/ATDHAR1/AOAT1/ANN1/ATPEPC1/GDPDL1/CLE12/AGT/ENO2/GOX1/ATPCB/AtBG2/CRK9/CORI3/PMDH2/BETA/UDG4/ATSBT4.12/LTP3/AtPRX71/ATBXL4/ANNAT2
                 Count
      GO:0009505    20
      GO:0048046    22
      
      GO基因集富集分析
      geneList = 2^e[,3]
      names(geneList)=as.character(ids[,3])
      geneList = sort(geneList, decreasing = TRUE)
      
      ego3 <- gseGO(geneList     = geneList,
                    OrgDb        = org.At.tair.db,
                    ont          = "CC",
                    nPerm        = 1000,
                    minGSSize    = 100,
                    maxGSSize    = 500,
                    pvalueCutoff = 0.05,
                    verbose      = FALSE)
      head(ego3)
      

      enrichGO和gseGO這兩個都用得比較多。

      可視化
      barplot(ggo, drop=TRUE, showCategory=12)
      dotplot(ego2)
      


      自己的分析略顯簡單,沒有過多的細化,后面隨著學習的深入再來補充。注:第2步提取gene_feature這一步沒有必要,因為軟件可以自動識別feature!2019.8.10


      reference

      Statistical analysis and visualization of functional profiles for genes and gene clusters:
      http://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html
      Analyzing RNA-seq data with DESeq2:
      http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

      因水平有限,有錯誤的地方,歡迎批評指正!

      bk
      下一篇:沒有了
    推土機2019:ImageCombiner - Java服務端圖片合成工具,好用! 萊布尼茨:【從零開始擼一個App】Fragment和導航中的使用 等不到的口琴:億級流量架構之資源隔離思路與方法 TOP生物信息:一文學會常規轉錄組分析 Twittytop:我的2020之路 我愛睡蓮:keepalived-1.3.5+MHA部署mysql集群 程序員養成日記:mysql一張表到底能存多少數據? 等你歸去來:java線程池趣味事:這不是線程池 Linyiwei:C++算法代碼――Tuna Linux中執行shell腳本的4種方法總結 關于ios配置微信config出現驗簽失敗的問題解決 ASP.NET Core擴展庫之Http通用擴展庫的使用詳解 Java在Excel中添加水印的實現(單一水印、平鋪水印) 基于UDP協議實現聊天系統 R語言中ifelse、which、%in%的用法詳解 constant,Java中定義常量(Constant)的幾種方法 float,CSS Float(浮動),CSS Float(浮動)用法 fixed,JavaScript fixed() 方法,JavaScript fixed() 實例 sprintf,sprintf 函數用法詳細注解 html代碼 php網站,php網站搭建步驟,PHP環境搭建教程 調試js,JS調試,js調試工具 xml是什么.xml格式文件,xml怎么用 jsswitch語句,JavaScript Switch 語句,JavaScript Switch用法 html5 教程,什么是 HTML5,HTML5都有什么功能 css,通過JS修改CSS樣式 網站建設哪,什么是網站建設?網站建設的常見要素有哪些? html網站,前端html網站的發布過程 document.cookie,使用document對象操作cookie javascript 數組,JavaScript數組去掉重復數據總結
    欧美日韩免费无码