- 博客(36)
- 收藏
- 关注
RSS订阅原创 RNA-seq上下游分析snakemake流程
使用fastp处理fastq文件，在使用START比对到基因组同时得到raw count，使用非冗余外显子长度作为基因的长度计算FPKM、TPM，同时也生成了CPM的结果。非冗余外显子长度计算可以参考之前的推文对定量结果质控使用生信技能树的三张图（PCA、树状图、热图）。使用python版的DEseq2对组间做差异分析（火山图和MA图）。
2024-06-12 21:41:10
512
原创 使用snakemake编写生信分析流程
后来才知道，reason不是推测的意思，而是名词原因的意思，这一步为什么会执行，因为输出文件不在指定的位置，换言之，如果我们跑完fastp_se后中断了snakemake流程，下次在接着跑流程，是不会跑fastp_se这一步的，因为这一步运行后输出了正确的文件results/trimmed/GSM6001951_L3.fastq。文件，虽然很长，其实就是一个判断你输入内容，然后交给fastp去执行的python脚本，所以我们需要按照作者的要求提供输入和输出文件名字，以及适当的额外参数。
2024-06-10 19:09:31
838
原创 Julia编程11：变量作用域 Scope of Variables
Julia有全局变量作用域和局部变量作用域，函数或者一些结构体、循环体如for等是否内部是局部环境可以参照下表。soft意思是只要不重名就行，hard必须是一个独立的局部变量。hard scope第一种情况很符合直觉，函数运行后x还是123，函数内的x默认是局部变量the hard scope rule isn't affected by anything in global scope第二种情况，使用global关键字声明x是全局变量，那么就修改了全局变量的x，因为全局变量中已经有了x，第二次赋值就
2024-06-04 22:57:27
586
原创 NMF（非负矩阵分解）分子分型
key_gene <- NMF::extractFeatures(result2, 0.5) # 提取关键基因。genelist ：基因名向量，基因名需要在exp的行名中，如果为空则使用全部基因。结果怎么看https://mubu.com/doc/C4gVcgp-G0。exp：表达矩阵，标准化过，但是不能有负数，行是基因，列是样本。method: 最常用的三种brunet、lee、snmf/r。ubuntu上，编译需要。n_run：运行次数。
2024-06-02 21:52:25
716
原创 Julia编程10：类型系统（Types）
Julia中给出了另一种方法，Vector{<:Real} 或者Vector{T} where T<:Real，这种类型实际上包括了多种类型，Vector{Int64}、Vector{Float64}、Vector{Int32}等等。基于类型系统的多重派发是Julia的核心，但是类型系统也挺复杂的，以下只是我的简单理解。typeof([1,2.0,3]) # 容器类型， Vector{Float64}Any是根类型，所有其他类型都是他的子类，如果函数的参数不指定类型则为Any类型。
2024-05-30 22:45:52
458
原创 Julia编程09：异常处理
println("√", xi, ": 平方根函数定义域异常")print("√", xi, ": 平方根函数其它异常")平方根函数，如果参数x<0，则报错"negative x not allowed"#利用短路逻辑，如果k<0就执行后边的语句 抛出错误。## √-2: 平方根函数定义域异常。无论如何最后都要执行的程序。## √a: 平方根函数其它异常。catch 异常类型变量名。# @assert相当于。
2024-05-28 20:45:34
501
原创 Julia编程08：控制流Control Flow
if 后边只能是布尔值即(true 或 false)，和Python不同可以放None表示False。如果i不是3的整倍数，那么就不运行println(i)，直接进入下一次循环。第一种相当于i和j的笛卡尔积，第二种是对应位置组合，如果位数不相同，那么只保留短的。a if test else b #如果test为真则返回a,否则返回b #注意R的ifelse是向量化的，也就是说test可以是多个布尔值。ifelse(test,a,b) # 向量化。嵌套for循环nested for。
2024-05-26 21:24:43
421
原创 Julia编程01：Julia语言介绍
参考https://zhuanlan.zhihu.com/p/585529899的代码，他所用的数字是n=40，但是我用julia算第40个斐波那契数只用了0.2s，感觉不太对，所以用了n=50，导致R和Python计算不出来，所以又折中算了n=45的情况。在2020上半年，因为疫情无法返校，所以在家待了半年，期间学习一点了R语言、Python、Julia、linux和C语言，只是学习基础语法并没有项目练习，因此返校半年后差不多都不记得了，现在重新捡起Julia丰富下当时写的笔记发到公众号来。
2024-05-24 20:25:21
343
原创 使用Arrow管理数据
在之前的数据挖掘：是时候更新一下TCGA的数据了推文中，保存TCGA的数据就是使用Arrow格式，因为占空间小，读写速度快，多语言支持（我主要使用的3种语言都支持）
2024-05-22 20:25:50
790
1
原创 文献阅读004:单细胞RNA测序揭示TME分析有助于化疗耐药膀胱癌患者治疗
scRNA-seq对别PDX靶向药（红色）和对照组（蓝色），发现处理组相较于对照组 HRAS表达水平下降（转录组水平和蛋白水平）；肿瘤细胞聚为2类，ClusterB更像是耐药细胞，但是对照组的一些细胞也分到了ClusterB中（蓝色），这说明耐药细胞一开始就存在于肿瘤样本中。肿瘤细胞分为3个cluster，其中cluster3的HRAS表达水平最低，IGFBP7, MDK, and B2M基因表达水平最高，休眠细胞最多。#1，#2，#3样本HE染色、WES和WTS测序。#3 样本：化疗耐药并且发生肺转移。
2024-05-20 20:11:31
358
原创 基因组实战03: WGS toy example
分析的数据是大肠杆菌，因为基因组小，适合拿来快速跑通整个流程。# WGS(全基因组测序)/WES(全外显子组测序)/WGRS(全基因组重测序)分析与可视化教程。ID：Read group identifier，必须唯一，一般为样本名；# 方法分享 | 全外显子测序分析的标准流程。#GATK4.0和全基因组数据分析实践。# GATK4全基因组数据分析最佳实践。# 重新合并过滤后的SNP和Indel。#第2节 FASTA和FASTQ。#第5节 理解并操作BAM文件。#第4节 构建WGS主流程。
2024-05-18 20:22:27
450
原创 Python基本统计分析
对于大样本，且频数表中每个单元格的期望频数都比较大（一般要求大于 5），可以不进行连续性校正。alternative可以选two-sided（默认，OR可能>1,也可能<1）, less（OR<1）, greater（OR>1）spss.fisher_exact这个函数的输入只能是2X2的二维列联表,R中的fisher.test输入可以不是2X2列联表。如果观察总例数 n 小于 40，或者频数表里的某个期望频数很小（小于 1），则需要使用 Fisher 精确概率检验。两个以上的变量交叉分类的频数分布表。
2024-05-14 19:42:35
497
原创 python获取网页表格数据
所以read_html是依靠lxml等库根据HTML语法找到表格位置，并转换为DataFrame。read_html 返回的是列表（a list of DataFrame）没有学过爬虫，好奇是read_html怎么做到的，怎么解析网页的。需要网页中的基因（Gene Symbol），一共371个。网页中的表格html语法大概如下。使用pandas读取网页表格。th: 定义表格的表头。tr: 定义表格的行。td: 定义表格单元。
2024-05-13 18:11:08
793
原创 所有物种基因Symbol别名转换为最新Symbol
在数据分析中会经常出现感兴趣的基因不在矩阵中，可能的原因是没有测到或旧版Symbol。因此需要找到旧版Symbol（Alias别名）和最新Symbol（Current Symbol）之间的对应关系。从NCBI下载最新的基因信息https://ftp.ncbi.nih.gov/gene/DATA/gene_info.gz。bq.tl.current_symbol可以把（表达）矩阵中的Symbol变为最新版。第二个参数Symbol与Alias对应关系文件路径。第一个参数数据框（index为Symbol）
2024-05-12 20:34:45
386
原创 python免费压缩PDF文件
去https://developer.ilovepdf.com/signup注册一个账户，然后把pulic key复制保存一下。每个月可以压缩250次PDF，对大多人是够的。如果extreme模式也不能压缩到需要的大小，可以编辑PDF，把矢量图换成像素图再进行压缩。'extreme'可以替换为'low', 'recommended'（默认）在终端设置pulic key（第一次需要设置，之后就不需要了）## 按装pdf-compressor包。Python环境准备。关于iLovePDF。
2024-05-11 18:35:35
227
原创 安装一个不存在的R包
需要用到DNA甲基化芯片的注释包IlluminaHumanMethylation450k.db，但是用BiocManager安装会报错“is not available”.搜到了生信技能树的文章《如果这个R包真的不存在了肿么办》，但只谈到包被删除了之前是存在的...修改IlluminaHumanMethylation450k.db/R/zzz.R中的源码。安装完才发现有替代的包 "FDb.InfiniumMethylation.hg19"...# 探针与基因名的对应关系。如何安装被移除的R包呢。
2024-05-10 20:23:15
1076
原创 数据挖掘：从表达谱芯片原始数据（CEL）到探针注释
range(exprs, na.rm = TRUE) # 1.889917 14.620563 不超过50不需要log2转化。CEL文件：探针的信号值和定位信息，是Affymetrix公司的芯片原始数据。# 得到探针对应的基因名字。3. limma标准化。4. 获取探针注释文件。0. 下载原始 数据。
2024-05-09 19:59:12
580
原创 数据挖掘：是时候更新一下TCGA的数据了
TCGA在去年更新之后提供了Count、TPM、FPKM三种格式的mRNA表达量数据，同时提供了ensembl gene ID、基因名、基因类型，因此有必要更新一下数据了。下载microRNA表达量数据。下载mRNA表达量数据。下载其他几种数据的函数。TCGA数据版本信息。
2024-05-08 20:51:18
1015
原创 DNA甲基化芯片分析02: DNA甲基化芯片基础知识基础知识
有的科学家觉得，DMR这样的区域还不够显著，DNA上的甲基化出现变化，可能是绵延几千位点的！而且只会在基因以外的区域，但是这些基因以外的区域发生变化，却会导致基因的表达发生变化。DMR：（连续的差异片段）一个连续不断都比较长的差异片段，科学家们觉得，这样的连续差异片段，对于基因的影响会更加明显，只找这样的片段，可以使得计算生物学的打击精度更为准确，也可以让最终找出来的结论数据更少，便于实险人员筛选。可以看到，位于编码基因上的探针最多，其次是位于基因间区的探针。可以看到，位于open sea的探针是最多的。
2024-05-07 19:36:11
500
原创 文献阅读003：单细胞肿瘤浸润性髓系细胞的异质性(综述)
综述主要介绍了Mast、Dendritic、Monocyte、Macrophage、Neutrophils细胞，其中Neutrophils不是很详细，因此只总结了前4中细胞的亚群marker和特征。幕布文档链接: https://www.mubucm.com/doc/7V2nk7N6B_1 密码: 9adx。表达精氨酸酶1、CD206、CD163、IL-4R、TGF-β1和血小板衍生生长因子(PDGF)促进Th2、抑制Th17、与CD4+ T互作促进肿瘤。促进血管生成、肿瘤生长和转移。
2024-05-06 21:24:32
1083
原创 使用极狐向github同步代码
由于不可抗力，github时断时续，而gitee有需要实名认证，所以作为替代品jihulab可以作为git同步代码，同时可以向gitlab和github同步代码。点击绿圈刷新便能同步，之后https://jihulab.com/BioQuest/scki更新后会自动向github推送更新，实现同步功能。进入https://jihulab.com/BioQuest/sckit在左侧设置中选择仓库。在https://github.com/settings/tokens生成token。
2024-05-05 17:31:47
244
原创 文献阅读02: 抑制CD39改变肿瘤微环境
流式结果表明CD39i增加CD27+ CD11b+和CD27-CD11b+ NK 细胞比例，但是减少 CD27-CD11b- and CD27 + CD11b- NK 细胞比例，这表明主要增加了成熟NK细胞；以上结果表明，CD39i 治疗增加NK细胞浸润，释放Xcl1从而招募cDC1，cDC1通过Il15, Il12b, Cd80 ，Cd86激活cycling CD8 + T细胞，同时也被免疫检查点抑制。图K说明在蛋白水平上tumor有更多T细胞，图LM说明CD39与CD8+ T细胞比例相关而与数量无关。
2024-05-04 18:37:56
910
原创 机器学习实战01:手动计算ROC-AUC
x轴是FPR，y轴是TPR，曲线上的每个点就对应着一组（FPR，TPR）坐标，所以我们的任务就是计算出所有的（FPR，TPR）坐标然后用线把他们连接起来就形成了ROC曲线，而AUC可以通过曲线下面积计算而来。计算AUC，AUC的定义是曲线下面积，按道理可以计算面积就行，但是如果样本较多，则会变成一条近似的曲线，计算了太大，因此有更好的方法计算AUC，比如。首先手动创建一组预测值和对应的真实值，0一般为负类，1为正类，而且正类多设置为研究中较为关心的标签，比如把1设置为肿瘤，或者疾病。FPR=0/5=0。
2024-05-03 19:49:52
857
原创 03.Python Dash网页开发：多页面网站制作
需求：写一个多网页的网站，包括header、footer、菜单包括主页home、博客blog（外部链接到博客）、about（自我介绍页面）每个页面都需要用dash.register_page(name, path='/') 声明是一个页面，path='/'设置路径为根路径。其中第一个APP是shiny APP，用外部网址链接到shiny网站，其他APP是DASH 示例APP。如果在子页面写callback需要把app.layout和@app.callback中的app去掉。
2024-05-02 19:11:54
925
原创 DNA甲基化芯片分析01: 使用methylumi和limma分析27K DNA甲基化芯片数据
27K的数据是很老的芯片数据，但是客户有需求就要找方法分析，主流的DNA甲基化芯片R包minfi和champ都只支持450K和850K的芯片。所以在bioconductor中搜索到了methylumi这个包，可以从idat读数据，经过质控得到beta值矩阵，之后用limma做差异分析。可以参考这篇文章[https://www.ncbi.nlm.nih.gov/pubmed/32375395]
2024-05-01 17:11:19
22
原创 文献阅读001: 由I类HLA转录缺失导致的联合免疫治疗的获得性癌症耐药性
背景Merkel cell carcinoma：梅克尔细胞癌又名皮肤小梁状癌、原发性皮肤神经内分泌癌、皮肤原发性小细胞癌及皮肤APUD瘤。HLA：MHC基因产物在不同细胞表面表达，通常称之为MHC分子或主要组织相容性抗原，对抗原递呈和免疫信号传递起关键作用。人的MHC又叫HLA。HLA抗原根据不同的基因位点分为一类、二类和三类抗原。一类抗原:HLA－A，B，C 位点二类抗原：HLA－DR，DQ，DP 位点三类抗原：补体摘要。
2024-04-30 19:52:38
696
原创 02:机器学习实战：最小二乘法
因此，比如我们的两个基因表达量相关性过高，就会造成多重共线性，建模中解决多重共线性的方法是可以对损失函数加上惩罚项，比如拉索回归、岭回归和弹性网络。可以看到，模型使用一条直线拟合散点，损失函数就是真实值到预测值的欧式距离的和，求损失函数最小值即是使直线尽可能拟合到更多的点。普通线性回归中的目标便是求的w和b两个参数，w其实是weight的简写，意为自变量的权重。：所给矩阵为行阶梯型矩阵，则矩阵中每一行的第一个不为零的元素的左边及其所在列以下全为零。是需要求得的未知数，而特征。此时的，w和x都是矩阵，，
2024-04-29 22:00:38
185
2
原创 02.Python Dash网页开发：网页有哪些元素组成与数据流
dcc.Input(id='my-input', value='王', type='text')中的value默认值为'王'，当在网页中输入新的字符后，value的值也会更新，更新后的value值传入update_output_div函数，返回的结果传递到my-output的children中，从而在网页中显示出来。如下所示，页面有两行，第一行有1列，宽度都是12，第二行分为3列，每列宽度是4（width=4）以下是一个简单的，交互式APP，我们输入字符，网站立刻返回 hello 字符。
2024-04-28 15:34:52
828
原创 01.Python Dash网页开发：环境配置和初试
我一直使用的是micromamba，因为比conda速度快，语法和conda一样，其中Dash网站所需要的4个包名字是dash开头，其他包是平时数据分析所需要用的；Dash类似R语言中的Shiny包，可以使用纯Python代码而不需要学习HTML、CSS、JavaScript语言就可以快速搭建一个网站，dash-bootstrap-components是Dash的拓展，提供了很多特性。这里使用的是dbc官网的案例，模仿Shiny包使用KMeans给iris数据集聚类。先不用管代码怎么写的，先跑起来。
2024-04-26 05:48:27
672
原创 教你如何优雅地管理R包
using函数是我写在.Rprofile中的函数，因此每次打开R就能使用，using的功能是一次加载多个包，并且使用了suppressPackageStartupMessages函数，因此不会显示加载包过程中的信息。使用pak管理R包，可以从Bioconductor、CRAN、Github、本地、URL安装R包，解决了R包安装需要多个不同R包去安装的问题。把代码框中的代码复制.Rprofile文件中，如果.Rprofile不存在则需要新建；写在.Rprofile文件中。使用using 加载R包。
2024-04-24 19:26:59
303
原创 单细胞测序 | marker基因列表怎么找？
单细胞转录组数据中比较核心且难以解决的事情是细胞注释，而细胞注释（Marker-based annotation）的质量又很大一部分取决于marker的数量和质量，我们常常苦于从何处寻找以及如何使用这些marker，上篇文章里提到的自动化注释工具ScType是将CellMarker和PanglaoDB数据库进行整合过滤后作为参考，对于我们来说难点在于如何实现其整合思想（关键在于各数据库中对细胞类型的命名缺乏统一和规范），以及如何支持更多的物种。[大笑R][满月R][微笑R]需要找我们 专研生信。
2024-04-22 22:36:23
529
空空如也
空空如也
TA创建的收藏夹 TA关注的收藏夹
TA关注的人