如何進行R-ChIP數據分析之準備篇

來自專欄生物科研

這應該是全網能搜索到的第一篇介紹R-ChIP數據分析的中文教程。

寫這篇文章的目的是想通過重複別人文章的分析提高下自己的分析水平。

本篇教程所用代碼更新在我的GitHub上，地址為https://github.com/xuzhougeng/R-ChIP-data-analysis

教程會對第一篇介紹R-ChIP技術文章"R-ChIP Using Inactive RNase H Reveals Dynamic Coupling of R-loops with Transcriptional Pausing at Gene Promoters"里的所有分析進行重複。

選擇這篇文章進行重複的理由有三點:

一：最近要探索R-loop數據分析流程

二：這篇文章的通訊作者是大牛，Xiang-Dong Fu

三：這篇文章將分析所用代碼都託管在https://github.com/Jia-Yu-Chen

背景知識

高通量數據分析其中在流程上大同小異，大多數會跑流程卻不知道如何分析的人，更多缺乏的是生物學基礎。在閱讀文獻後，我整理下幾個我覺得和後續數據分析有關的幾個知識點:

• R-loop是一種RNA/DNA三鏈結構體，與基因組穩定性和轉錄調控有關。
• 通過電鏡觀察，R-loop大小在150~500bp之間。
• 硫酸氫鹽測序(bisulfate sequencing)表明R-loop主要出現在基因啟動子的下游。
• R-loop所在非模板鏈(又稱編碼鏈)具有很強的序列偏好性，計算方式為(G-C)/(G+C)

R-loop的高通量分析方法目前都是依賴於S9.6抗體捕獲RNA/DNA雜合體，然後超聲打斷或酶切，如果後續對DNA進行測序，那就是DRIP-seq(DNA:RNA immunoprecipitation [DRIP] sequencing)，如果後續對RNA逆轉成的cDNA繼續測序，那就是 [DRIPc]-seq(DNA:RNA immunoprecipitation followed by cDNA conversion)。 然而酶切的解析度不夠，超聲又容易破壞脆弱的R-loop結構，於是就導致目前很多文獻報道有矛盾。

這篇文章就開發了一種新方法，基於RNase H的體內R-loop譜檢測策略。作者構建一種沒有催化活性的RNASE H1(原始的RNase H會同時靶向和基因組和線粒體基因組，因此在N端加上了NLS用於定位，在C端加上V5標籤用於IP ），RNASEH1與RNA/DNA結合，超聲打碎，用anti-V5抗體進行染色體免疫共沉澱(ChIP)。隨後RNA/DNA雜合體轉換成雙鏈DNA(ds-DNA), 之後便是鏈特異性測序。

關於鏈特異性測序，推薦拜讀鏈特異性測序那點事

「和普通的RNAseq不同，鏈特異性測序可以保留最初產生RNA的方向」

文章的分析內容

後續的分析部分會嘗試完成文章中出現的分析。

證明R-ChIP方法的效率和可靠度

• 找到D210N特異而WKKD非特異的基因
• peak的正負鏈比對情況

序列偏好性和基因組分布

• 比較narrow Peak和broad peak的peak大小（boxplot)
• broad peak是否包含narrow peak（韋恩圖），且一方特異的peak的信號值也弱
• G/C skew， 統計距離summit的ATCG比例，統計背景和R-loop相比，G二聚體、G-三聚體、G四聚體..的出現比例
• G-rich 在非模板鏈出現次樹顯著增加
• peak在基因組上的分布(promoter proximal regions, +- Kb from TSS). 比較不同基因不同部分的差異

與S9.6的R-loop結果系統性比較(基於K562 cells)

• 比較DRIP-seq和R-ChIP的peak大小
• 看兩者的peak的overlap
• 在JUN座位上的R-ChIP, DRIP-seq, DNase, H3K4me1, H3K4me2,H3K4me3, H3K27ac

準備分析環境

軟體部分

文章中"Software and Algorithms"這部分列出了分析主要所用的軟體，加上下載SRA數據所需工具和一些常用軟體，一共要安裝的軟體如下:

• SRA Toolkit: 數據下載工具
• Bowtie2: 比對工具
• SAMtools: SAM格式處理工具
• BEDtools: BED格式處理工具
• MACS2: 比對後找peak
• R: 統計作圖
• Ngsplot: 可視化工具
• Deeptools: BAM文件分析工具, 可作圖。

軟體安裝部分此處不介紹，畢竟如果你連軟體安裝都有困難，那你應該需要先學點Linux基礎，或者去看生信必修課之軟體安裝

分析項目搭建

使用 mkdir創建項目文件夾，用於存放後續分析的所用到的數據、中間文件和結果

1. mkdir -p r-chip/{analysis/0-raw-data,index,scripts,results}

個人習慣，在項目根目錄下創建了四個文件夾

• analysis: 存放原始數據、中間文件
• index: 存放比對軟體索引
• scripts: 存放分析中用到的腳本
• results: 存放可用於放在文章中的結果

後續所有的操作都默認在 r-chip下進行，除非特別說明。

數據下載

根據文章提供的GEO編號(GEO: GSE97072)在NCBI上檢索, 按照如下步驟獲取該編號下所有數據的元信息, 我將其重命名為"download_table.txt"然後上傳到伺服器 。

微信不支持過大的GIF圖，想看GIF去我的簡書

然後在 scripts下新建一個腳本 01.sra_download.sh

1. #!/bin/bash
2. sra_files=${1?missing input file}
3. sra_ids=$(egrep -o SRR[0-9]{5,9} $sra_files)
4. mkdir -p analysis/log
5. for id in $sra_ids
6. do
7. prefetch $id >> analysis/log/download.log
8. done

在項目根目錄下運行 bash scritps/01.sra_download.sh ownload_table.txt

下載的數據默認情況下存放在 ~/ncbi/public/sra, 需要用 fastq-dump解壓縮到 analysis/0-raw-data. fastq-dump的使用說明見Fastq-dump: 一個神奇的軟體

新建一個腳本，叫做 02.sra_uncompress.sh，存放在 scripts文件下，代碼如下

1. #!/bin/bash
2. file_in=${1?missing input file}
3. SRR_IDS=$( egrep -o SRR[0-9]{5,9} ${file_in})
4. OUT_DIR="analysis/0-raw-data"
5. mkdir -p ${OUT_DIR}
6. for id in $SRR_IDS
7. do
8. echo "fastq-dump --gzip --split-3 --defline-qual + --defline-seq @$ac-$si/$ri $id -O ${OUT_DIR} &" | bash
9. done

然後用 bash02.sra_uncompress.sh download_table.txt運行。

如果空間不夠用，可以用 egrep-oSRR[0-9]{5,9}download_table.txt|xargs-i rm~/ncbi/public/sra/{}.sra刪除SRR數據。

注意：這是單端測序，所以每個SRR只會解壓縮出一個文件

此外還需要下載human genome (hg19)的bowtie2索引，用於後續bowtie2比對。

1. curl -s ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/hg19.zip -o index/hg19.zip &
2. cd index
3. unzip hg19.zip

從索引提取每條染色體的長度

1. bowtie2-inspect -s index/hg19 | cut -f 2,3 | grep ^chr > index/hg19.chrom.sizes

來源：生信媛

ps:如有侵權，請聯繫小編進行刪除。