学术咖

Mindblown: a blog about philosophy.

  • preprocessCore分位数标准化

    library(“preprocessCore”) library(“ggplot2”) data <- read.table(“mouse_gc.xls”,sep=”\t”,header = TRUE) rownames(data) <- data$GeneID #rownames(data) <- data[,1] #data_mat <- data.matrix(data[,-1]) quantile_normalisation <- function(df){ df_rank <- apply(df,2,rank,ties.method=”min”) df_sorted <- data.frame(apply(df, 2, sort)) df_mean <- apply(df_sorted, 1, mean) index_to_mean <- function(my_index, my_mean){ return(my_mean[my_index]) } df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean) rownames(df_final) <- rownames(df) return(df_final) } new_data <- quantile_normalisation(data[,2:19]) boxplot(data) boxplot(new_data) write.table(new_data,…

  • 跨物种比较转录组

    近几年几篇高分文献对跨物种比较转录组方法的应用及研究领域，感兴趣的可自行查看原文献： Cell_2019_跨物种单细胞测序揭示小胶质细胞基因表达谱进化特征 Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱 Science_2018_跨人类和猕猴大脑发育的时空转录组差异 Cell_2021_African lungfish genome sheds lighton the vertebrate water-to-landtransition Cell_2021_人类大脑进化的发育调控机制 Cell_2021_原始辐鳍鱼类基因组的解析发现登陆的遗传基础已经在硬骨鱼类祖先出现 看完上述文献可知和种内比较转录组相比，跨物种（种间）比较转录组要复杂的多，主要体现在背景基因总数、测序深度、管家基因表达谱和可能的其他因素差异等；分析流程主要分为两步：1、直系同源基因的查找。此处大部分文献使用 one-to-one orthologs （及双向最优比对），如果物种跨度较大可选较为宽松的阈值，期望值1e-10，覆盖度50%；如果物种亲缘关系较为接近就选较严的阈值。也有文献先获得直系同源基因列表，一个gene list及为一个基因（matagene），该基因的表达量为gene list所有基因表达量的和（基因拷贝数差异）。2、合并后基因表达矩阵的标准化。标准化方法各异，方法之一：首先在种内计算出FPKM、RPKM、TPM等，种内TMM标准化，然后根据背景基因集合并数据得到新的表达矩阵再使用TMM标准化。也有直接quantile normalization using the R package preprocessCore（具体R实现可参考下篇文章），标准化方法的选择对结果的影响很大，可自行尝试适合自己的标准化方法。另外，不同数据来源需要去除批次效应，不然进行PCA和样本相关性分析时会将不同数据来源的样本聚到一块；数据的适当缩放也是必要的，一些聚类方法会将极值带来的影响放大化，使结果不准确。下图是Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱分析流程。

  • 数据相关性分析

    简介：评价两组数据之间的相关性，有皮尔森（pearson）相关系数，斯皮尔曼（spearman）相关系数和肯德尔（kendall）相关系数。在这三大相关系数中，spearman和kendall属于等级相关系数亦称为“秩相关系数”，是反映等级相关程度的统计分析指标。相关系数的绝对值越大，相关性越强，相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。 Pearson（皮尔逊）相关系数：皮尔逊相关也称为积差相关（或积矩相关）是英国统计学家皮尔逊于20世纪提出的一种计算直线相关的方法。适用于： (1)、两个变量之间是线性关系，都是连续数据。 (2)、两个变量的总体是正态分布，或接近正态的单峰分布。 (3)、两个变量的观测值是成对的，每对观测值之间相互独立。 Spearman Rank（斯皮尔曼等级）相关系数：在统计学中，斯皮尔曼等级相关系数以Charles Spearman命名，并经常用希腊字母ρ（rho）表示其值。又称秩相关系数，是利用两变量的秩次大小作线性相关分析，对原始变量的分布不作要求，属于非参数统计方法，适用范围要广些。斯皮尔曼等级相关是根据等级资料研究两个变量间相关关系的方法。它是依据两列成对等级的各对等级数之差来进行计算的，所以又称为“等级差数法”。斯皮尔曼等级相关对数据条件的要求没有积差相关系数严格，只要两个变量的观测值是成对的等级评定资料，或者是由连续变量观测资料转化得到的等级资料，不论两个变量的总体分布形态、样本容量的大小如何，都可以用斯皮尔曼等级相关来进行研究。对于服从Pearson 相关系数的数据亦可计算Spearman 相关系数，但统计效能要低一些。Pearson 相关系数的计算公式可以完全套用 Spearman 相关系数计算公式，但公式中的x 和y 用相应的秩次代替即可。 Kendall Rank（肯德尔等级）相关系数：在统计学中，肯德尔相关系数是以Maurice Kendall命名的，并经常用希腊字母τ（tau）表示其值。用于反映分类变量相关性的指标，适用于两个分类变量均为有序分类的情况。对相关的有序变量进行非参数相关检验； 取值范围在-1-1 之间，此检验适合于正方形表格；肯德尔(Kendall)W 系数又称和谐系数，是表示多列等级变量相关程度的一种方法。适用这种方法的数据资料一般是采用等级评定的方法收集的，即让K 个 评委（被试）评定N 件事物，或1 个评委（被试）先后K 次评定N 件事物。等级评定法每个评价者对N 件事物排出一个等级顺序，最小的等级序数为1 ，最大的为N ，若并列等级时，则平分共同应该占据的等级，如，平时所说的两个并列第一名，他们应该占据1 ，2 名，所以它们的等级应是1.5, 又如一个第一 名，两个并列第二名，三个并列第三名，则它们对应的等级应该是1,2.5,2.5,5,5,5, 这里2.5 是2,3 的平均，5 是4,5,6 的平均。 R包计算相关性系数： library(ggplot2) library(ggpubr) data <- read.table(“111.xls”,sep=”\t”,header = TRUE)#列是特征（两个特征），行是物种名 ggplot(data=data, aes(x=genome, y=Tes))+geom_point(color=”red”)+stat_smooth(method=”lm”)+stat_cor(data=data, method = “pearson”)+theme_classic() #stat_cor(data=dat, method…

  • MUMmer使用及后续绘图

    MUMmer被广泛用于大片段序列的比对，如染色体共线性分析。 nucmer [options] delta-filter -i 80 -l 1000 -r -q out.delta > out.rq.delta #比对率大于80%，对比长度大于1000 #-r: 仅保留每个reference在query上的最佳位置,允许多条reference在query上重叠 #-q: 仅保留每个query在reference上的最佳位置,允许多条query在reference上重叠 show-coords out.rq.delta > out.coords grep -P “zfCh04\s+” out.coords|awk ‘{print $12,$13}’ |sort |uniq -c #查看reference单个染色体zfCh04和query不同染色体的匹配区域的数量 grep -P “zfCh04\s+gcCh04” out.coords|awk ‘{print $12,$1,$2,$13,$4,$5}’ > zfCh04_gcCh04.txt sed -i ‘s/ /\t/g’ zfCh04_gcCh04.txt #R作图 使用RIdeogram包，可参考RIdeogram：染色体数据可视化的R包 install.packages(‘RIdeogram’) require(RIdeogram) cc <- read.table(“111.xls”,sep=”\t”,header = TRUE,stringsAsFactors =…

  • ATAC-seq分析流程（三）

    前边介绍了使用DiffBind进行ATAC-seq数据的差异分析，DiffBind进行差异分析时调用了edgeR或者DESeq2包；相对于DiffBind流程我更习惯普通转录组分析流程，因为两者本质上是一样的，区别在于ATAC-seq要根据结果的peak文件自己制备gtf文件，peak对应的gtf文件和gene是不一样的，不同样本间也可能存在差异。 cat *.narrowPeak > merged.peaks sort -k1,1 -k2,2n merged.peaks > merged.peaks.sorted bedtools merge -i merged.peaks.sorted > bedtools.merged.sorted.bed awk -v var=ATAC_seq ‘{print $1″\t”var”\texon\t”$2″\t”$3″\t.\t+\t.\tgene_id \”peak_”NR”\”; transcript_id \”peak_”NR”\”;”}’ bedtools.merged.sorted.bed > bedtools.merged.atac.gtf featureCounts -T 18 -p -t exon -g gene_id -a bedtools.merged.atac.gtf -o ATAC_counts.txt *last.bam sed 1d ATAC_counts.txt | cut -f 1,7-15 > rawATAC_counts.matrix 一个peak相当于一个gene，得到的原始矩阵可以使用DESeq2包分析，筛选出共有特有或者差异的peaks等进行peaks的注释。

  • ATAC-seq分析流程（二）

    六、Deeptools可视化 #deeptools安装conda install deeptools gtf2bed < ~/gc_genome/17.Geta/1.geta/out.gtf > genome.bed for i in `cat ../samples.txt`; do echo “bamCoverage -p 5 –normalizeUsing RPKM -b $i.last.bam -o $i.last.bw”; done > bam2bw.list ParaFly -c bam2bw.list -CPU 9 computeMatrix reference-point –referencePoint TSS -p 15 -b 10000 -a 10000 -R ../genome.bed -S ../*bw –skipZeros -o matrix1_test_TSS.gz –outFileSortedRegions regions1_test_genes.bed plotHeatmap -m matrix1_test_TSS.gz…

  • BS-seq分析流程（一）

    短序列的质控都可以使用trimmomatic，这里不多做介绍，得到的clean data可做下面分析 一、比对和甲基化位点提取（Bismark） Bismark安装及使用 git clone https://github.com/FelixKrueger/Bismark.git #conda install -c bioconda bismark #Genome Preparation /opt/biosoft/Bismark-0.23.0/bismark_genome_preparation –parallel 40 –verbose ./ #bismark alignment for i in `cat ../samples.txt` do echo “/opt/biosoft/Bismark-0.23.0/bismark –parallel 40 –genome ./ –phred33-quals -1 ../$i.1.fastq -2 ../$i.2.fastq” done > bismark.list ParaFly -c bismark.list -CPU 4 #deduplicate for i in `cat ../samples.txt` do echo “/opt/biosoft/Bismark-0.23.0/deduplicate_bismark…

  • ATAC-seq分析流程（一）

    首先根据今年发表在Genome Biology上一篇综述总览分析流程： 作者建议的预处理分析流程为：FastQC-trimmomatic-BWA-MEM-ATACseqQC。 一、使用Fastqc对fastq 文件的测序质量统计 for i in `ls *fastq`; do echo “fastqc -t 4 -o ./ $i”; done > fastqc.sh nohup ParaFly -c fastqc.sh -CPU 18 & 二、Trimmomatic质控（可使用 cutadapt,AdapterRemoval v2,Skewer 和 trimmomatic 等软件） The forward and reverse adapters are slightly different. We will also trim low quality bases at the ends of the reads…

  • Upset plot：进阶版Venn plot实现多元素可视化

    if (!requireNamespace(“BiocManager”, quietly = TRUE)) install.packages(“BiocManager”) BiocManager::install(“UpSetR”) library(UpSetR) pathway=read.table(“carp_combined_id1.xls”,header=T,sep=”\t”) upset(pathway,nsets = 9, nintersects = 50,mb.ratio = c(0.5, 0.5), order.by = “freq”, decreasing = c(TRUE,FALSE), queries = list(list(query=intersects, params=list(“datra”, “caaur”,”meamb”,”onmac”,”ctide”,”cycar”,”sigra”,”sians”,”sirhi”), color=”red”, active=T))) #详细参数查看官网 #行列转置awk ‘{for(i=0;++i

  • R画折线图

    library(ggplot2) library(gcookbook) data

Got any book recommendations?

Get In Touch