R-loop数据分析之R-ChIP(样本间BAM比较和可视化)

最新推荐文章于 2024-06-22 09:37:20 发布
最新推荐文章于 2024-06-22 09:37:20 发布
阅读量635 收藏 2
点赞数 1
文章标签： r语言
原文链接：https://yq.aliyun.com/articles/650485
版权
样本间相关性评估

上一步得到各个样本的BAM文件之后，就可以在全基因组范围上看看这几个样本之间是否有差异。也就是先将基因组分成N个区间，然后用统计每个区间上比对上的read数。

脚本scripts/07.genome_bin_read_coverage.sh如下

#!/bin/bash

set -e
set -u
set -o pipefail

samples=${1?missing sample file}
chromsize=${2:-index/hg19.chrom.sizes}
size=${3:-3000}
ALIGN_DIR="analysis/2-read-align"
COV_DIR="analysis/3-genome-coverage"

mkdir -p ${COV_DIR}

exec 0< $samples

while read sample
do
    bedtools makewindows -g $chromsize -w $size | \
      bedtools intersect -b ${ALIGN_DIR}/${sample}.flt.bam -a - -c -bed > ${COV_DIR}/${sample}.ReadsCoverage
done

准备一个输入文件存放待处理样本的前缀，然后运行脚本bash scripts/07.genome_bin_read_coverage.sh samples_rep.txt

HKE293-D210N-Input-Rep1
HKE293-D210N-Input-Rep2
HKE293-D210N-Input-Rep3
HKE293-D210N-V5ChIP-Rep1
HKE293-D210N-V5ChIP-Rep2
HKE293-D210N-V5ChIP-Rep3

最后将得到的文件导入到R语言中进行作图，使用的是基础绘图系统的光滑散点图(smoothScatter)。

files_list <- list.files("r-chip/analysis/3-genome-coverage","ReadsCoverage")
files_path <- file.path("r-chip/analysis/3-genome-coverage",files_list)

input_rep1 <- read.table(files_path[1], sep='\t')
input_rep2 <- read.table(files_path[2], sep='\t')
input_rep3 <- read.table(files_path[3], sep='\t')
chip_rep1 <- read.table(files_path[4], sep='\t')
chip_rep2 <- read.table(files_path[5], sep='\t')
chip_rep3 <- read.table(files_path[6], sep='\t')

pw_plot <- function(x, y, 
                    xlab="x",
                    ylab="y", ...){
  log2x <- log2(x)
  log2y <- log2(y)
  smoothScatter(log2x,log2y,
                cex=1.2,
                xlim=c(0,12),ylim=c(0,12),
                xlab=xlab,
                ylab=ylab)
  text(3,10,paste("R = ",round(cor(x,y),2),sep=""))
}

par(mfrow=c(2,3))
pw_plot(chip_rep1[,4], chip_rep2[,4],
        xlab = "Rep 1 (Log2 Tag Counts)",
        ylab = "Rep 2 (Log2 Tag Counts)")

pw_plot(chip_rep1[,4], chip_rep3[,4],
        xlab = "Rep 1 (Log2 Tag Counts)",
        ylab = "Rep 2 (Log2 Tag Counts)")

pw_plot(chip_rep2[,4], chip_rep3[,4],
        xlab = "Rep 1 (Log2 Tag Counts)",
        ylab = "Rep 2 (Log2 Tag Counts)")

pw_plot(input_rep1[,4], input_rep2[,4],
        xlab = "Rep 1 (Log2 Tag Counts)",
        ylab = "Rep 2 (Log2 Tag Counts)")

pw_plot(input_rep1[,4], input_rep3[,4],
        xlab = "Rep 1 (Log2 Tag Counts)",
        ylab = "Rep 3 (Log2 Tag Counts)")

pw_plot(input_rep2[,4], input_rep3[,4],
        xlab = "Rep 2 (Log2 Tag Counts)",
        ylab = "Rep 3 (Log2 Tag Counts)")
img_c61437c8226f0eabf60a169416d27c73.jpe
correlationship

这种多个BAM文件之间相关性衡量，其实也可以用deepTools的plotCorrelation画出来，但是我觉得应该没有R语言画 的好看。

BigWig可视化

由于同一个样本间的BAM文件具有很强的相关性，因此可以将这些样本合并起来转换成bigwig格式，这样子在基因组浏览器（例如IGV, UCSC Browser, JBrowse）上方便展示。

如下的代码的目的就是先合并BAM，然后转换成BigWig，拆分成正链和反链进行保存

#!/bin/bash

set -e
set -u
set -o pipefail

samples=${1?please provied sample file}
threads=${2-8}
bs=${3-50}

ALIGN_DIR="analysis/2-read-align"
BW_DIR="analysis/4-normliazed-bw"

mkdir -p $BW_DIR

exec 0< $samples
cd $ALIGN_DIR

while read sample
do
    bams=$(ls ${sample}*flt.bam | tr '\n' ' ')
    if [ ! -f ../$(basename ${BW_DIR})/${sample}.tmp.bam ]; then
    echo "samtools merge -f -@ ${threads} ../$(basename ${BW_DIR})/${sample}.tmp.bam $bams  &&
          samtools index ../$(basename ${BW_DIR})/${sample}.tmp.bam" | bash
    fi
done

cd ../../
exec 0< $samples

cd ${BW_DIR}
while read sample
do
    bamCoverage -b ${sample}.tmp.bam -o ${sample}_fwd.bw -of bigwig  \
      --filterRNAstrand forward --binSize ${bs} --normalizeUsing CPM --effectiveGenomeSize 2864785220 \
      --extendReads 150 -p ${threads} 2> ../log/${sample}_fwd.log
    bamCoverage -b ${sample}.tmp.bam -o ${sample}_rvs.bw -of bigwig  \
      --filterRNAstrand reverse --binSize ${bs} --normalizeUsing CPM --effectiveGenomeSize 2864785220 \
      --extendReads 150 -p ${threads} 2> ../log/${sample}_rvs.log
    rm -f ${sample}.tmp.bam ${sample}.tmp.bam.bai
done

得到的BW文件你可以在IGV上初步看看，比如说检查下文章Figure 1(E)提到的CIRH1A基因

img_d63d0e2b348f70e6bae8ea933e41516d.jpe
CIRH1A

发表文章时肯定不能用上图，我们可以用R的Gviz进行展示下

library(Gviz)

#下面将scale等track写入tracklist
tracklist<-list()
itrack <- IdeogramTrack(genome = "hg19", chromosome = 'chr16',outline=T)
tracklist[["itrack"]]<-itrack

# 读取BigWig
bw_file_path <- "C:/Users/DELL/Desktop/FigureYa/R-ChIP/"
bw_file_names <- list.files(bw_file_path, "*.bw")
bw_files <- file.path("C:/Users/DELL/Desktop/FigureYa/R-ChIP/",
                      bw_file_names)

tracklist[['D210-fwd']] <- DataTrack(range = bw_files[3],
                              genome="hg19",
                              type="histogram",
                              name='D210 + ',
                              ylim=c(0,4),
                              col.histogram="#2167a4",
                              fill.histogram="#2167a4")
tracklist[['D210-rvs']] <- DataTrack(range = bw_files[4],
                                     genome="hg19",
                                     type="histogram",
                                     name='D210 - ',
                                     ylim=c(4,0),
                                     col.histogram="#eb1700",
                                     fill.histogram="#eb1700")

tracklist[['WKDD-fwd']] <- DataTrack(range = bw_files[9],
                                     genome="hg19",
                                     type="histogram",
                                     name='WKDD + ',
                                     ylim=c(0,4),
                                     ylab=2,
                                     col.histogram="#2167a4",
                                     fill.histogram="#2167a4")
tracklist[['WKDD-rvs']] <- DataTrack(range = bw_files[10],
                                     genome="hg19",
                                     type="histogram",
                                     name=' WKDD -',
                                     ylim=c(4,0),
                                     col.histogram="#eb1700",
                                     fill.histogram="#eb1700",
                                     showAxis=TRUE)

#写入比例尺
scalebar <- GenomeAxisTrack(scale=0.25,
                            col="black",
                            fontcolor="black",
                            name="Scale",
                            labelPos="above",showTitle=F)
tracklist[["scalebar"]]<-scalebar

# 画图
plotTracks(trackList = tracklist,
           chromosome = "chr16",
           from =  69141913, to= 69205033,
           background.panel = "#f6f6f6",
           background.title = "white",
           col.title="black",col.axis="black",
           rot.title=0,cex.title=0.5,margin=10,title.width=1.75,
           cex.axis=1
           )
img_38ebfc28a44cd43226f7dd1ef79cf731.jpe
Gviz

后续用AI修改下坐标轴，就几乎和原图差不多了。此外可能还要调整之前脚本的bin size， 使得整体更加平滑

看图说话： 由于WKKD在D210N的基础上继续突变了几个位点，使得原本只丧失了RNASEH1的催化活性的D210N进一步丧失了结合到RNA/DNA杂合体的能力，在实验中就可以作为 负对照。上图就是其中一个有代表性的区间。

weixin_33766168
关注
  • 1
    点赞
  • 踩
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
2020.9.17丨Chip-seq结果可视化之peak检测（下）
09-17 3634
这一部分是使用deeptools对样品进行相关性分析以及主成分分析，同时从peak中去挖掘motif。我使用的工具是MEME-ChIP，MEME是一个工具系列，挖掘motif的工具比较丰富，MEME、DREME、TomTom、MEME-ChIP，其中MEME-ChIP可以同时调用其他几个工具进行综合分析，比较方便。 deeptools是一个很好用的深度分析工具，中文版使用手册可以让你快速上手（虽然翻译有些直，但竟然看得懂！）。在进行相关性分析和主成分分析之前，需要对样品数据进行一个综合统计，deeptoo
vizbam:BAM 可视化
05-30
已弃用 这个软件已经过时了。 威兹巴姆 显示SAM / BAM对齐的 Java 库。 项目状态：测试版 作者：皮埃尔·林登鲍姆博士。 @yokofakun 依赖关系 Java 1.7 皮卡德 阿帕奇蚂蚁 截图 显示像samtools tview这样的 BAM（项目一部分） 命令行工具 可视化命令行 用法：VizBamCmdLine [选项] 像 samtools tview 一样打印 SAM 对齐 选项 描述 格式=字符串 输出格式 默认值：文本。 可以将此选项设置为 'null' 以清除默认值。 参考=文件 输入索引引用 必需。 输出=文件 输出文件 默认值：null。 输入=文件 要处理的 BAM 文件。 必需的。 位置=字符串 要处理的位置。 语法是 "chrom:position" 染色体必须存在于参考中。 默认值：空。 INTERVALFILE =文件 包含间隔
bam 可视化：samtools tview 详细解释
qq_28723681的博客
09-27 9422
bam可视化软件大家可能熟悉的是IGV，然而，IGV对于大多数linux用户来说并不友好，而samtools tview可以很好满足该需求。 话不多说直接上命令行： samtools tview -p chr1:3128088 NA12878.bam hg38.fasta -p 指定染色体的位置，tview从指定的位置开始显示 NA12878.bam 比对结果bam文件,需要构建索引（NA12878.bam.bai） hg38.fasta 比对时使用的fasta文件，如果不提供则第一行会显示为N. 需要
推荐项目：trackplot——在R中快速可视化bigWig文件的利器
最新发布
gitblog_00056的博客
06-22 315
推荐项目：trackplot——在R中快速可视化bigWig文件的利器 项目地址:https://gitcode.com/PoisonAlien/trackplot 项目简介 面对基因组数据分析中的大体积数据，trackplot是一个响应迅速且易于上手的R脚本工具，专注于生成类似IGV风格的轨迹图（也称为位点图）、轮廓图和热图，专门针对bigWig文件。它简化了复杂的数据可视化过程，为生物信息学研...
零基础小白笔记8 | 将bam文件转换为bw文件并进行可视化
tmrtmr___的博客
01-02 3297
2.使用工具：deeptools是一个用于深度分析测序数据的工具集，其中包含多个子命令，其中bamCoverage可以将BAM 文件转换成 BigWig 格式的输出文件。1.bw文件，即bigwig文件，是一种常用的用于存储基因组测序数据的文件格式，它允许高效存储大量基因组范围的信号数据，并支持快速的数据检索与可视化。（2.2）设置参考基因组为mm10（本次实验动物模型），选择染色体及区域，最左侧轴设置显示区域的大小；所以bw文件在转换完成后需要先下载到本地才能在IGV上进行可视化分析；
ChIPseeker包对ChIP-seq peaks可视化
qq_27390023的博客
06-19 1484
ChIPseeker软件包实现了以下功能：检索峰值附近最近的基因、注释峰值的基因组区域、估计芯片峰值数据集之间重叠重要性的统计方法，以及整合地理数据库，以便用户将自己的数据集与存储在数据库中的数据集进行比较。比较可以用来推断合作规则，因此可以用来生成假设。实现了几个可视化功能，以总结峰实验的覆盖范围、结合到TSS区域的峰的平均轮廓和热图、基因组注释、到TSS的距离以及峰或基因的重叠。...
DROMPAplus：ChIP-seq管道工具，用于对多个ChIP-seq样本进行质量检查，标准化，统计分析和可视化
02-26
DROMPA（DRaw和观察多种浓缩曲线和注释）是一种ChIP-seq管道工具，可满足各种需求，包括质量检查，分析和多个ChIP样品的可视化。 DROMPAplus用C ++编写，具有许多有价值的功能。 DROMPAplus： 接受多种地图文件...
ChIP-seq数据分析.pdf
12-23
ChIP-seq 数据分析 ChIP-seq（Chromatin Immunoprecipitation sequencing）是一种高通量的表观遗传学实验技术，用于研究蛋白质与 DNA 的相互作用。该技术通过结合特异性抗体，捕获与蛋白质结合的 DNA 片段，然后...
hictoolsr:R软件包，其中包含用于Hi-C数据分析的有用工具的集合
03-26
`hictoolsr` 是一个专为Hi-C数据分析设计的R软件包，它提供了丰富的工具集，旨在帮助生物信息学家和研究人员处理、分析和可视化Hi-C实验产生的大量数据。Hi-C技术是一种用于研究细胞核内染色质相互作用的高通量方法...
NanoR:R程序包用于分析和比较纳米Kong数据-开源
04-29
NanoR是一款专为R编程语言设计的开源软件包，它主要致力于分析和比较一维MinION和GridION X5测序数据，特别是纳米孔（Nanopore）测序技术所产生的Kong数据。Nanopore测序技术是近年来发展迅速的新型基因组测序方法，...
Python-NGS数据分析工具代码
08-10
6. **可视化**：为了便于理解和解释结果，数据通常需要进行可视化，Python的`matplotlib`、`seaborn`、`plotly`等库提供了强大的绘图功能。 "ngstools"可能包含了上述部分或全部功能，提供了Python接口来简化NGS...
transrate-tools:读取bam文件并汇总有关序列的读取映射信息
05-06
转换工具 用于处理bam文件的命令行工具。 用C ++编写并使用Bamtools。 建造 transrate-tools使用C ++ 11功能，因此您至少需要安装g ++ 4.7。 在OSX上，您可以使用brew install gcc49安装最新的gcc。 请注意，在OSX上，您始终需要使用选项-DCMAKE_CXX_COMPILER=$(which g++-4.9)来告诉cmake在哪里可以找到gcc。 确保使用子模块克隆： $ git clone --recursive git@github.com:cboursnell/transrate-tools.git 并且您将需要安装cmake。 接下来，构建bamtools： 在Linux上： $ cd bamtools $ mkdir build $ cd build $ cmake .. $ make $ cd ../..
deepTools对ChIP-seq数据可视化
qq_27390023的博客
06-25 1247
deepTools是一套python工具，特别为高效分析高通量测序数据而开发，如ChIP-seq、RNA-seq或MNase-seq。 参考：https://deeptools.readthedocs.io/en/develop/content/tools/computeMatrix.html http://deeptools.readthedocs.org/en/latest/content/example_gallery.html#normalized-chip-seq-signals-and-p
ChIP-seq分析流程(基于linux系统)
生信小博士的博客
12-16 2285
1）自己实验获得的原始测序文件（.fastq.gz），测序文件格式为.fastq格式，而.gz格式是压缩格式，为了节省存储空间。 可以将测序文件修改为自己样品的名字 2）从NCBI的GEO数据库中下载sra文件，需要转换为.fastq文件。 $ wget -c ‘下载地址’ #wget下载数据，-c参数为断点续传 从GEO数据库中下载的.sra文件 #如果是单端测序文件就只需要fastq-dump命令将sra文件转换为fastq测序文件。 $ fastq-dump SRR10485695.sra #如果
CHIP-seq流程学习笔记（10）-使用intersect()函数提取ChIP-seq和RNA-seq中的共有基因
垚垚爸的博客
11-07 738
参考文章：R中常用函数使用说明 # 利用intersect函数对RNA_seq和ChIP_seq中的交集gene_symbol进行提取。 # 参考文章：https://blog.csdn.net/woodcorpse/article/details/80494605?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160471856319195264746932%2522%252C%2522scm%2522%253A%252220140713.
chip画图
qq_39306047的博客
06-07 1606
library(ChIPseeker) library(ggplot2) library(dplyr) library(org.Mm.eg.db) require(TxDb.Mmusculus.UCSC.mm10.knownGene) txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene#TxDb文件 #f <- getSampleFiles()[[4]]#.bed文件,这里用包里的测试数据 DM1 <- "C:\\Users\\admin\\Desk.
2020.9.15丨Chip-seq结果可视化之peak检测（上）
09-15 3632
macs2运行参数 macs2 callpeak -t K1_ChIPed_S1_L007_R1.bam -c K1_Input_S5_L007_R1.bam -f BAM -g mm -n K1 -B -q 0.01 -t -c 实验组和对照组结果 ​-f 输入文件格式 -g 参考基因组有效大小，人类选择hs，也可以根据基因组大小直接输入数值 -n 输出前缀 -B 输出bdg格式文件，可以上传到UCSC生成峰图 -q q值，默认0.05 -p p值，未校正值 导入到R中
samtools从fastq到bam再到bigwig（bw）
Gentlezzx的博客
12-03 4285
首先需要使用bowtie2和samtools将fastq序列比对到参考基因组上，生成.bam文件： 单末端： "bowtie2 -p 10 -x mm10-U input.fq | samtools sort -O bam -@ 10 -o - > output.bam 双末端： bowtie2 -p 10 -x mm10-1 input_1.fq -2 input_2.fq | samtools sort -O bam -@ 10 -o - > output.bam 这里的mm10是bow
对分子模拟轨迹数据的分析绘图
基博的博客
09-02 4471
主要是用gnuplot和origin绘制分子模拟分析常用的图
RNA-seq生信分析R语言代码
09-27
RNA-seq的生信分析可以使用R语言进行。下面是RNA-seq分析的大体步骤和对应的R语言代码： 1. 数据质控： - 质量评估和过滤：使用软件如FastQC和Trimmomatic对测序数据进行质量评估和过滤。 ```R library("ShortRead") library("Rsubread") # 质量评估 qc <- FastqFile("sample.fastq") quality <- qa(qc) # 过滤 filtered <- filterFASTQ(qc, output_file="filtered.fastq", trimVector=TRUE) ``` 2. 数据比对： - 将测序数据比对到参考基因组或转录组上，使用软件如HISAT2和STAR。 ```R library("Rsubread") # 比对到参考基因组 index <- buildindex("genome.fasta", "genome_index") aligned <- align(index, "filtered.fastq", output_file="aligned.bam") ``` 3. 确定基因表达水平： - 使用软件如featureCounts和HTSeq将比对结果转换为基因表达矩阵。 ```R library("Rsubread") # 生成基因表达矩阵 counts <- featureCounts("aligned.bam", annot.ext="gene.gtf", isPairedEnd=FALSE) ``` 4. 差异基因表达分析： - 使用软件如DESeq2和edgeR对基因表达矩阵进行差异基因分析。 ```R library("DESeq2") # 差异分析 dds <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info, design=~condition) dds <- DESeq(dds) results <- results(dds) ``` 以上是RNA-seq生信分析的大体步骤和对应的R语言代码。请根据实际情况对代码进行相应的调整和优化。如果需要更详细的代码和解释，请参考相关的生物信息学教程和文献。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值