博文

R 正则表达式

 jg <- read_table('~/Documents/sjbd.csv',col_names = F) %>%    mutate(X1=str_replace(X1,"^\"","1")) %>%    mutate(X1=str_replace(X1,",\\\\N","")) %>%    filter(!str_detect(X1,"^\\,")) %>%    filter(!str_detect(X1,"^\\(")) %>%    separate(X1,c('A','B','C','D','E','F','G','H','I','M','K','N','O'),sep=",") jg <- read.table('~/Documents/sjbd.csv') 

解决 Windows11 下 wsl 2使用 matplotlib.pyplot 画图没反应的问题

图片
  获取wls2主机地址 并 添加 .bashrc文件中  Dynamically export the  DISPLAY   environment variable in WSL2; # add to ~/.bashrc, and source it. export DISPLAY=$(cat /etc/resolv.conf | grep nameserver | awk '{print $2}'):0 2. 下载并打开  VcXsrv  on Windows. It’s  import  to check the “Disable access contr ol” to allow connection from WSL2; 3. pycharm中设置 Setting the DISPLAY environment variable in PyCharm;

python链接oracle数据库

import pandas as pd import numpy as np import janitor import oracledb # 链接方式一 params = oracledb.ConnectParams ( host = "192.168.30.48" , port = 1521 , service_name = "JZDB1" ) con = oracledb.connect ( user = "ipvsdb" , password = "######" , params =params ) # 链接方式二 dsn= f' { "ipvsdb" } / { "#####" } @ { 1521 } :192.168.30.48/JZDB1' con=oracledb.connect ( dsn ) sqldq = ''' select dzbm,dzmc from sys_xzqh_zzjg ''' xzqh = pd.read_sql ( sqldq, con ) con.close ()

解决PyCharm中Python Console的中文乱码问题

  进入Settings->Build, Execution, Deploymer->Console->Python Console Environment variables中添加PYTHONIOENCODING=UTF-8 Settings->Editor->File Encoding 均选择UTF-8

RNA-seq 流程

 #cuffdiff pipeline FASTQ-reads -> mapping(HISAT2,tophat2) ->BAM ->cuffdiff ->DEGs(FPKM) #DESeq2 pipepline FASTQ-reads -> mapping(HISAT2,tophat2) ->BAM ->Gene Raw Reads Count ->find DEGs in R (No FPKM)

RNA-seq脚本备份

 hisat2-build建立索引 hisat2-build ./data/Saccharomyces_cerevisiae.R64-1-1.dna_rm.toplevel.fa ./data/yeast_ref hisat2将read比对到参参考基因组由于sam格式的文本文件过大,一般将其转换为bam格式的文件 hisat2 -x ./data/yeast_ref -U ./data/SRR1916152.fastq | samtools view -bS -| samtools sort - -o ./data/EV_3.bam hisat2 -x ./data/yeast_ref -U ./data/SRR1916153.fastq | samtools view -bS -| samtools sort - -o ./data/EV_4.bam hisat2 -x ./data/yeast_ref -U ./data/SRR1916154.fastq | samtools view -bS -| samtools sort - -o ./data/DNMT3B_2.bam hisat2 -x ./data/yeast_ref -U ./data/SRR1916155.fastq | samtools view -bS -| samtools sort - -o ./data/DNMT3B_3.bam hisat2 -x ./data/yeast_ref -U ./data/SRR1916156.fastq | samtools view -bS -| samtools sort - -o ./data/DNMT3B_4.bam 建立bam文件的索引 samtools index ./data/EV_3.bam samtools index ./data/EV_4.bam samtools index ./data/DNMT3B_2.bam samtools index ./data/DNMT3B_3.bam samtools index ./data/DNMT3B_4.bam 安装HTSeq软件,计算每个基因中reads的数目 htseq-count -f bam -r pos ./data/EV_3.bam ./data/Saccharomyces...

计算fpkm

setwd("E:\\read_count") library(tidyverse) library(magrittr) library(hablar) library(DESeq2) library(GenomicFeatures)  ##计算fpkm #根据参考基因组获得每条reads的exon txdb <- makeTxDbFromGFF("E:\\read_count\\Saccharomyces_cerevisiae.R64-1-1.106.gtf",format="gtf") exons_gene <- exonsBy(txdb, by = "gene") exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))}) n=t(as.data.frame(exons_gene_lens)) #传入Count矩阵和基因长度 FPKM <- function(data_count,gene_len) {   lens <- nrow(data_count)   colnames <- data_count %>% colnames()   for (i in 1:ncol(data_count)) {     mapped_reads <- sum(data_count[1:lens,i])#计算每个样本的mapped reads数,由于data_count文本最后几行是综述数据,需要舍去,因此选择58721行之前的基因数据。     FPKM <- data_count[1:lens,i]/(10^-9*mapped_reads*gene_len[1:lens,1])#计算FPKM值,10^-9是单位转换后的换算,不能放在分母外的原因是mapped_reads*n[1:58721,1]的乘积很大会报错,因此需要先缩小倍数再进行计算     data_count = data.frame(data_count[1:lens,],FPKM)#添加到data_count表格后面i   }   data_coun...