生信小tips

CHarLoTTe♚

已于 2024-07-10 11:47:42 修改

阅读量2.5k

点赞数 3

分类专栏： RNA-seq chipseq 文章标签： r语言

于 2021-12-30 09:55:50 首次发布

本文链接：https://blog.csdn.net/Yuan_CHarLoTTe/article/details/119489222

版权

chipseq 同时被 2 个专栏收录

6 篇文章 0 订阅

订阅专栏

RNA-seq

5 篇文章 2 订阅

订阅专栏

说明：做生信过程中有很多小细节需要注意，留个备忘录方便后面查找命令

对染色体命名列内容修改

在这里插入图片描述
第一列为染色体号，如果要修改成chr格式

 sed -i -e 's/^/chr/g' （）.narrowPeak

在这里插入图片描述
第一列为chr，删掉chr改为数字

sed -i 's/^chr//g' 1_summits.bed

在这里插入图片描述

将sra文件转成fastq格式

yuanye@DNA:~/projects/dongfeng/ChIP_seq/2021_7_29/data/fastq$ ~/software/sratoolkit.2.11.0-ubuntu64/bin/fastq-dump --split-3 SRR…….1

将fastq格式压缩

yuanye@DNA:~/projects/dongfeng/ChIP_seq/2021_7_29/data/fastq$ gzip -c SRR…….1_1.fastq > SRR…….1_1.fastq.gz

将fq.gz解压

gunzip 0h-******_1.clean.fq.gz

多文件解压和质控

 gunzip *.fq.gz

fastqc -t 16 -o /Data/wu_lab/yuanye/projects/RNAseq/daizhongye/2021_8_11/QC/  *.clean.fq

提取某文件的两列

awk '{print $1,"\t"$2}' *.isoforms.results > gene_trans.txt

** 一定是isoforms文件

去除ENSEMBLID 中的小数点

DESeq_Txi4 <- read.csv("DESeq-Txi42.csv",header = T)
head(DESeq_Txi4)
DESeq_Txi4$gene_id=unlist(str_split(DESeq_Txi4$gene_id,"[.]",simplify=T))[,1]

构建condition

condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))

将bam文件转为bigwig文件（用STAR比对后）###

出现ERRO:LI……s1.0.0得时候
conda install samtools=1.13
即可解决

然后对bam文件进行sort，index

多文件数据取交集合并（R）

setwd("F:/rnaseq/dongfeng/2021_8_16/data2/heatmap")
library(dplyr)
A <- read.csv("Acountid.csv",header = F)
B <- read.csv("Bcountid.csv",header = F)
C <- read.csv("Ccountid.csv",header = F)
OVER <- read.csv("overlapAB_BC.csv",header = F)
OVERA <- left_join(OVER, A, by ='V1')
OVERAB <- left_join(OVERA, B, by ='V1')
OVERABC <- left_join(OVERAB, C, by ='V1')

结果
在这里插入图片描述

绘制heatmap

注：去掉重复名字的数据

install.packages("pheatmap")
library(gplots)
library(pheatmap)
heatABC <- read.csv("OVERABC.csv",header = T)
data <- read.csv("OVERABC.csv", header = T, row.names = 1)
pheatmap(as.matrix(data), scale='row', border_color=NA, width=9, height=50, fontsize=9, fontsize_row=10,cluster_col=F, filename="heatmap2.pdf" )
pheatmap(data, scale='row', border_color=NA, width=9, height=50, fontsize=9, fontsize_row=10,cluster_col=F, filename="heatmap3.pdf" )

R语言学习笔记（3）——left_join() right_join() inner_join() full_join()

homer注释时motif.bed文件的染色体命名需要以“chr”开头

环境MOTIF

findMotifsGenome.pl motif.bed hg19 /Data/wu_lab/sunxiaoyu/environment/R8_motif/ -len 8,10,12

添加chr
cat file |awk -F "\t" '{OFS="\t"} {if (1==1)1="Chr"1}1' >file_chr
或
awk '{print "chr"$0}' xxx.sorted.gtf > xxx.sorted.gtf

去掉chr

sed -i 's/chr//' xxx.sorted.gtf

将转换完的ID与差异表达基因表合并

https://zhuanlan.zhihu.com/p/358168557#:~:text=GSEA （Gene Set EnrichmentAnalysis）

RNA比对软件选择

由于分析的 RNA ，那么就存在可变剪切，对于存在可变剪切的 mapping 用 Tophat 或者 Tophat 的升级工具 HISAT2 更合适

√查看后台运行程序

cat nohup.out

√ H3K36ME2找broad peak

√ 查看conda安装过的软件

在这里插入图片描述

√ Linuxubuntu64版本

conda 环境激活失败

wu_lab@bio:~$ cd /Data/wu_lab/anaconda3/bin/
 chmod 777 activate
 source ./activate

就会出现base环境

如果要默认进入base环境则

cd /Data/wu_lab/miniconda3/bin/
chmod 777 activate
source ./activate
cd  ~
source ./.bashrc
source ./.profile

网页做GO和KEGG

kobas

http://kobas.cbi.pku.edu.cn/genelist/

书写csv文件出现以下错误

在这里插入图片描述文件内容为
保存整个文件会报错，但是我改为只保存counts内容，则解决了。

write.csv(y_T98G_GAS41sh4_1_vs_scr$counts,file="normolize_y_T98G_GAS41sh4_1_vs_scr.csv")

如果写xls文件
则

library(xlsx)
write.xlsx(xxx, "xx.xls",sheetName = "xxx",append = TRUE)

批量解压

for gz in *.gz; do gunzip $gz; done# 解压多个.gz文件
for tar in *.tar.gz;  do tar xvf $tar; done# 解压多个.tar.gz文件

Linux查看磁盘剩余空间大小

Linux系统下查看硬盘大小的方法：首先在Linux系统中打开终端；然后在终端中输入“df -h”命令即可查看硬盘大小，获取硬盘被占用了多少空间，目前还剩下多少空间等信息。 df命令可以获取硬盘被占用了多少空间，目前还剩下多少空间等信息，它也可以显示所有文件系统对i节点和磁盘块的使用情况。

https://www.jianshu.com/p/5eb30a71c2c0

详情

df [-ahikHTm] [目录或文件名]

选项与参数：

-a ：列出所有的文件系统，包括系统特有的 /proc 等文件系统；
-k ：以 KBytes 的容量显示各文件系统；
-m ：以 MBytes 的容量显示各文件系统；
-h ：以人们较易阅读的 GBytes, MBytes, KBytes 等格式自行显示；
-H ：以 M=1000K 取代 M=1024K 的进位方式；
-T ：显示文件系统类型, 连同该 partition 的 filesystem 名称 (例如 ext3) 也列出；
-i ：不用硬盘容量，而以 inode 的数量来显示

选项与参数：

-a ：列出所有的文件与目录容量，因为默认仅统计目录底下的文件量而已。
-h ：以人们较易读的容量格式 (G/M) 显示；
-s ：列出总量而已，而不列出每个各别的目录占用容量；
-S ：不包括子目录下的总计，与 -s 有点差别。
-k ：以 KBytes 列出容量显示；
-m ：以 MBytes 列出容量显示；

参考：Linux下内存查看命令
在这里插入图片描述

修改Linux进入的登录默认目录路径

要修改登录后默认目录，可以用vim编辑器，打开/etc/passwd ，找到相应的用户，修改倒数第一个冒号前面的目录即可，如下图所示。
在这里插入图片描述

修改文件后缀

例：将xls改为bed格式

mv DIPG17-vor_H3K27me3_loss.{xls,bed}

用xls转换出来的bed自动换列，如果用csv出来的没有换列

快捷命令

快捷方式：

Ctrl + U   剪切光标位置到行首的字符

Ctrl + E   回到行尾

Ctrl + A   回到行首

Ctrl + W   剪切一个单词

Ctrl + Y   粘贴命令行剪切的内容

Ctrl + K   剪切光标位置到行首的字符



一些命令：

cd -   # 返回上次的工作目录

ls -t   # 以时间排序列出当前目录下文件

ls -R   # 递归目录列出文件

ls -d   # 显示目录本身，而非目录下文件

ls ./*txt   # 列出当前目录下以txt结尾的文件

ls ../   # 列出上层目录的文件

ls -l   # 列出当前目录下文件的详细信息

mkdir -p   # 递归创建目录

mv   # 移动或重命名

cp -r   # 复制

rm -rfi   # 递归删除/不显示警告讯息/删除前询问用户

head/tail   # 查看文件前/后n行，默认10行。head常结合管道符用于控制输出行数。（管道符：|）

more   # 逐页查看，按空格翻页，按回车换行，q退出

less -S   # 单行显示

less -N   # 显示行号

zless   # 查看压缩文件

cat   # 查看文本文件内容，输出到屏幕 作者：一只小蛮要 https://www.bilibili.com/read/cv15760548?from=search&spm_id_from=333.337.0.0 出处：bilibili

「生信技能树」2021公益课(linux基础 & conda)

如果没有spike in 如何进行normalize

scale factor = 总reads数 / 测序深度

如果仅仅转bigwig则可利用RPKM/FPKM
RPKM的诞生是针对SE（单端）测序，FPKM则是在PE（双端）测序上对RPKM的校正。

比对后PCR去重

软件：samtools rmdup / picard

如何写bed文件

写bed文件可以先写成txt文件然后将文件后缀名改成bed文件即可

call peak 工具

1.macs2
2.sicer
3.epic2
4.HMMRATAC （可用于ATAC- seq call peak ）
5.使用MACS2执行H3K27ac、H3K4me1和H3K4me3信号的峰值调用。SICER用于调用H3K27me3数据的峰值。所有峰值呼叫的FDR阈值均为0.01。

搜索conda安装软件版本

conda search 软件名

批量修改文件名中的某字符

Linux rename命令

 rename "s/_val_1//" *

去除文件中全部带_val_1的字符
或者
找到所有以fq.gz后缀结尾的文件，将他们的“-”全部修改成“_”，但是一次只能修改一个，也就是h3-k27-me3，要修改两次

find  *.fq.gz | rename 's/-/_/'

查看后台运行命令以及终止

jobs进行查看
运行的命令有数字号码
然后 kill %号码即可终止（不需要加[ ]）
在这里插入图片描述

基因组名称以及大小

GRCm38=mm10；2652783500
GRCm39=mm39；
GRCh37=hg19； 2864785220
GRCh38=hg38。 2913022398

bash 快捷键们

在这里插入图片描述

mv移动文件

在这里插入图片描述

删除命令

在这里插入图片描述

vim光标移动

在这里插入图片描述

linux知识

Linux Tools Quick Tutorial
鸟哥的 Linux 私房菜 – 基础学习篇目录
 Welcome to Linux World

broad or narrow

在这里插入图片描述

bed 参考基因组文件下载

bed文件的下载：

hg19:https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/

mm10:https://sourceforge.net/projects/rseqc/files/BED/Mouse_Mus_musculus/

两个都选择RefSeq

【转录组入门】5：序列比对

log计算器

莫得感情计算器

某列全部+1

由于利用bedtools coverage出来的counts会出现某段有0count得情况，对后续处理数据产生影响。因此在counts列每个数都+1.

cat EZH2_WT_K119ub_position_EZH2_wt_K119ub.xls|awk '{print $4+1}' > EZH2_WT_K119ub_position_EZH2_wt_K119ub+1.xls

prefetch下载srr文件

prefetch --option-file data.txt

将txt里的逗号替换成换行符

sed ":a;N;s/,/\n/g;ta" a.txt > PRCa.csv

提取基因GTF文件染色体位置和基因名

用的是gencode 的gtf文件

awk '{if((NR>5) && ($2=="HAVANA") && ($3=="gene")) print $0}' gencode.v38.annotation.gtf | perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' > hg38geneposition.gft.bed

参考链接：生信编程11.根据gtf格式的基因注释文件提取人所有基因的染色体坐标

R软件包遇到非零退出怎么办？

解决R包“had non-zero exit status”安装报错

install.packages("stringr",type="binary")

R语言中交集，并集，补集，差集的方法汇总

如何获取sam文件中序列比对到哪条链，FLAG信息解读

在这里插入图片描述

 pigz -k *.fastq -p 20

-k 保留原始文件
-p 线程数

awk 匹配指定列包含的指定字符串并输出该行

目的:我要输出E2 所在行
在这里插入图片描述

 awk -F "\t" '$4 == "E2" {print}' ESC_50_chrhmm_segments.bed > E2.bed

E2在哪列$就几

bedtools coverage计算注意事项

bedtools coverage 计算时不要以0为起始
如： chr1 0 12345 否则下面的chr2/chr3就不会计算了即遇到0就会停止计算
有时候分完bin后可能会存在下一个染色体起始成为上一个染色体的终止，导致出现这种情况 chr2 9999 100 这种也是要去除的，否则不能计算counts

统计文件行数

 wc -l xxx.bed

将文件按行数平均分成几等份

linux下的split 命令（将一个大文件根据行数平均分成若干个小文件）

split -l 3253  BAP1_AID_WT_0h_GenesFPKMlow2_mm10StartEnd_postion_COUNTSOF_H2AZ.bed -d -a 4 BAP1_AID_WT_0h_GenesFPKMlow2_mm10StartEnd_postion_COUNTSOF_H2AZ

本次文件共13009行，想要平均分成4等分，所以每个文件内是约3253行

出现rlang包版本错误怎么办

R语言：解决报错"载入了名字空间‘rlang’ ..，但需要的是>= .."
首先在Rstudio中卸载相应的包（DESeq2.rlang），

remove.packages('DESeq2')

然后关掉Rstudio
在R中（不是Rstudio）
在这里插入图片描述

以及DESeq2的包
重新在Rstudio中加载即可

R中某列数乘一个数或除以一个数

如何仅将 R 数据框中的一列与数字相乘？

countsummary$Reads <- countsummary$Reads/2
countsummary$Percentage <- countsummary$Percentage*2

黑名单

Boyle-Lab/Blacklist

Linux：查找指定字符串，并将结果输出到指定文件

grep ‘1805’ CloudPayment.log | grep ‘1905’ > out.log

subset 进行数据筛选（&且 |或）

从零开始学R数据分析，数据筛选与提取

R中删除某列

在这里插入图片描述

conda 安装软件出现下列错误即解决方法

(base) wu_lab@bio:~$ conda update conda
Collecting package metadata (current_repodata.json): / WARNING conda.models.version:get_matcher(546): Using .* with relational operator is superfluous and deprecated and will be removed in a future version of conda. Your spec was 1.7.1.*, but conda is ignoring the .* and treating it as 1.7.1                                                  done
Solving environment: | ^Z
[7]+  已停止               ( "$CONDA_EXE" $_CE_M $_CE_CONDA "$@" )
(base) wu_lab@bio:~$ conda config --remove-key channels
(base) wu_lab@bio:~$ conda update conda
Collecting package metadata (current_repodata.json): done
Solving environment: done

在这里插入图片描述

tar文件解压

tar -xf xxxxxx.tar

下载NCBI参考基因组

1.先搜索想要的物种名称点击ASSembly
在这里插入图片描述
2.下载fna(fa)文件以及gtf文件

3.解压

tar -xf Ecoli_assemblies_genome_gtf.tar
tar -xf Ecoli_assemblies_genome_fasta.tar

cd到miniconda bin中的bowtie-build路径

 bowtie2-build /Data/wu_lab/reference/NCBI_Ecoli/ncbi-genomes-2023-10-10/GCF_000008865.2_ASM886v2_genomic.fna /Data/wu_lab/reference/NCBI_Ecoli/bowtie_index/

提取某个文件夹下的所有或指定后缀的文件名

 ls -1 | grep ".gz$" > 1.txt

linux将某目录下所有文件名写到一个文件里

log2FC R计算

#data=log2(data[,1:6]+1) #对基因表达量数据处理
 
#data <- as.matrix(data) #转变为matrix格式矩阵
 
#head(data

分组

在这里插入图片描述

df <- read.csv("percentage.csv")
df=data.frame(x=df$log2FCWT,y=df$log2FCKO)
for(i in 1:length(df$x))
{if (df$x[i]>=df$y[i]) df$z[i]="up" 
else df$z[i]="down"}
print(df)