跳转至

生物信息学

  • 分类于 生物信息学
  • 需要 1 分钟阅读时间

计算panel bed编码区长度

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#!/bin/bash
#@File    :   run.sh
#@Time    :   2023/11/21 09:32:17
#@Author  :   biolxy
#@Version :   1.0
#@Contact :   biolxy@aliyun.com
#@Desc    :   None

inputbed=$1

# export PATH=/data/biogonco/lixy/bedtools/bedtools2/bin:$PATH  # v2.30.1 不行
export PATH=/data/bioinfo_project/bioinfo_miniconda3/bin:$PATH  # v2.30.0 不行
# /data/biogonco/lixy/bedtools/bedtools2/bin/bedtools

bedtools merge -i ${inputbed} > merge.bed

bedtools summary -i merge.bed -g /mnt/nas_001/project/oncodemo/genome_and_annotation/hg19/gatk/b37/chrom.sizes > merge.bed.summary  # 统计的是 bed 区间的长度


# 去 gencode 下载 gencode.v41.annotation.gff3
# wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/GRCh37_mapping/gencode.v44lift37.annotation.gff3.gz
# /mnt/nas_001/project/oncodemo/genome_and_annotation/genome/gffread-0.12.7.Linux_x86_64/gffread gencode.v44lift37.annotation.gff3 -T -o gencode.v44lift37.annotation.bed
# awk '$3 == "CDS"' gencode.v44lift37.annotation.bed | sed 's/^chr//g' > coding_CDS.bed

# coding_CDS.bed 为编码区bed

bedtools intersect -a merge.bed -b /mnt/nas_001/project/oncodemo/genome_and_annotation/genome/coding_CDS.bed | sort -k1,1 -k2,2n |  bedtools merge -i - > panel.cds.bed

# panel.cds.bed 为 pane 的 编码区 bed
bedtools summary -i panel.cds.bed -g /mnt/nas_001/project/oncodemo/genome_and_annotation/hg19/gatk/b37/chrom.sizes > panel.cds.summary

注意

要注意bedtools的版本问题，我测试的 v2.30.0可以，v2.30.1 不行

  • 分类于 生物信息学
  • 需要 1 分钟阅读时间

准确度和灵敏度和特异性

定义

Precision 准确度：指找出来的突变中的真阳性有多少

Recall 召回率：指总突变集合有多少可以被找出来

公式

准确度（Precision）：

准确度是指在所有被预测为阳性（即突变）的样本中，实际为阳性（真阳性）的比例。换句话说，它衡量的是预测结果中真阳性的比例。其计算公式为：

\(\text{Precision} = \frac{TP}{TP + FP}\)

其中，TP是真阳性的数量，FP是假阳性的数量。

召回率（Recall）或 灵敏度/敏感度（Sensitivity）：

召回率是指在所有实际为阳性的样本中，被正确预测为阳性（即真阳性）的比例。它衡量的是测试方法捕捉到所有实际阳性样本的能力。其计算公式为：

\(\text{Recall} = \frac{TP}{TP + FN}\)

其中，TP是真阳性的数量，FN是假阴性的数量。

这两个指标通常用于评估分类模型的性能，特别是在医学和生物学研究中，它们帮助研究者理解一个测试或分析方法在识别阳性样本方面的准确性和完整性。在体细胞变异检测的背景下，这两个指标对于评估变异检测算法的性能至关重要。

特异性（Specificity）：

特异性是指在实际没有某种疾病的人群中，诊断测试能够正确排除非患者的能力。它衡量了测试对非患者的"特异性"，即测试能够准确地排除非患者的能力。 特异性的计算公式为：

特异性 = 真阴性（True Negative）/（真阴性 + 假阳性（False Positive））

\(\text{Specificity}= \frac{TN}{FP + TN}\)

注意： 这里不使用 特异性 这一概念，原因是在NGS数据calling中，难以确定真阴性位点的数量，特异性无法计算。

F1_score

F1_score，是评估分类模型性能的一种指标，特别是在二分类问题中。它结合了精确度（Precision）和召回率（Recall）两个指标来提供一个单一的评分，以衡量模型的整体性能。 F1_score 是精确度和召回率的调和平均数，公式为：

\(F1 = 2 \times \frac{\text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}\)

F1_score 的值范围是0到1，1表示完美的精确度和召回率，0表示模型性能很差。

在处理不平衡的数据集时，F1_score 特别有用，因为它同时考虑了精确度和召回率，而不是仅仅关注其中一个。这使得 F1_score 成为一个在多种情况下都相对平衡的性能度量标准。

  • 分类于 生物信息学
  • 需要 1 分钟阅读时间

生物信息学习技术路线(生信指月录)

1. Linux

推荐 Linux就该这么学 （0-5章节，其他的可选择性学习）

  1. 要求：
    1. 掌握 shell 基础语法
      1. Linux基础命令
      2. cd, ls, mkdir, pwd, time, df, cp, rm, sed, awk, wc, head, tail, more, history等
    2. 目录如下：
      1. 第0章 咱们先来谈谈为什么要学习Linux系统
      2. 第1章 动手部署一台Linux操作系统
      3. 第2章 新手必须掌握的Linux命令
      4. 第3章 管道符、重定向与环境变量
      5. 第4章 Vim编辑器与Shell命令脚本
      6. 第5章 用户身份与文件权限

2. Python

  1. 掌握 基础语法，条件判断，数据类型，简单的数据结构
  2. 熟悉常用模块：sys, os, pandas, json, Bio, numpy, seaborn, pysam, argparse等
  3. 使用Python 读写文件（xls, xlsx, txt等文件）
  4. 了解 迭代器，生成器，装饰器等概念，并描述其适用的场景
  5. python2 和 python3 中 range() 函数的区别
  6. 可以被next()函数调用并不断返回下一个值的对象称为迭代器：Iterator (ref: https://www.liaoxuefeng.com/wiki/1016959663602400/1017323698112640)
  7. 在Python中，这种一边循环一边计算的机制，称为生成器：generator (ref: https://www.liaoxuefeng.com/wiki/1016959663602400/1017318207388128) 熟悉 类 和 实例的概念，类函数，静态函数
  8. 可选：
    1. Python 多线程操作
      1. 使用 pymysql 操作数据库
      2. 使用 python-docx等操作数据库
      3. 常见的设计模式及其使用场景
      4. Python 数据结构与算法
      5. 陈斌老师 《数据结构与算法》B站课程
      6. 课件： 链接: https://pan.baidu.com/s/1srvCWOLsxEn3mOq1_TCmhQ 提取码: wtji

3. 生物信息学课程

  1. 了解相关概念, 专有名词
  2. 了解相关概念背后的 生物信息算法，统计原理
  3. 【 山东大学】生物信息学 B站课程
  4. 【北京大学】生物信息学：导论与方法 B站课程
  5. 统计学 ⅰ. https://www.yuque.com/biolxy/bioinfo/zgl1ml

4. BioStar 实操练习 https://www.biostarhandbook.com/

  • 可以看生信媛翻译的 BIOSTAR课程

5. 肿瘤外显子数据分析指南

  • https://www.yuque.com/biotrainee/wes

sad

  • 分类于 生物信息学
  • 需要 2 分钟阅读时间

获取启动子序列

获取启动子序列

  1. 获取物种基因组序列genome.fa文件
  2. 获取物种注释信息gff文件

输入参数： positionFile, genomeFa, outfa （pos坐标文件，包含四列 chrom start end strand）

使用Python脚本

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
'''
@File    :   Promoter.py
@Time    :   2021/08/15 20:16:07
@Author  :   biolxy 
@Version :   1.0
@Contact :   biolxy@aliyun.com
@Desc    :   None
@Usage   :   python Promoter.py pos.Lee.txt /home/lixy/soybean/other/glyma.Lee.gnm1.BXNC.genome_main.fna out.Lee.fa
'''
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO
import sys


class Promoter(object):
    @staticmethod
    def get_promoter_by_genome(positionFile, genomeFa, outfa, length=2000):
        """positionFile
        Chr20   35557820    35562522    +

        length: 2000
        """
        dict1 = {}
        pos_list = Promoter.get_pos_list(positionFile)
        for seq_record in SeqIO.parse(genomeFa, "fasta"):
            seqid = str(seq_record.id)
            for tmp_l in pos_list:
                if seqid == tmp_l[0]:
                    seq = seq_record.seq
                    s = int(tmp_l[1]) - length
                    e = int(tmp_l[1])

                    if '-' == tmp_l[3]:
                        s = int(tmp_l[2])
                        e = int(tmp_l[2]) + length
                    # print(s, e)
                    seq = seq[s:e]
                    if '-' == tmp_l[3]:
                        seq = seq.reverse_complement() # 反向互补
                    n_seqid = seqid + "_" + str(tmp_l[1]) + ":" + str(tmp_l[2]) + "_" + str(tmp_l[3])
                    dict1[n_seqid] = seq

        records = []
        for seqid in dict1:
            seq = dict1[seqid]
            rec1 = SeqRecord(seq, id=seqid, description="")
            records.append(rec1)
        SeqIO.write(records, outfa, "fasta")


    @staticmethod
    def get_pos_list(infile):
        list_ = []
        with open(infile, 'r') as ff:
            for line in ff:
                line = line.strip()
                tag = line.split("\t") # chrom start end strand
                list_.append(tag)
        return list_


if __name__ == "__main__":
    positionFile, genomeFa, outfa = sys.argv[1], sys.argv[2], sys.argv[3]
    Promoter.get_promoter_by_genome(positionFile, genomeFa, outfa)

  • 分类于 生物信息学
  • 需要 3 分钟阅读时间

由 reads count 数计算 TPM 和FPKM

由 reads count 数计算 TPM 和FPKM
  • 下载gff文件
  • http://venanciogroup.uenf.br/cgi-bin/gmax_atlas/download_by_pub.cgi 去这里下载reads count 数据
S
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
library("GenomicFeatures")
library("dplyr")
path1 = setwd("C:/Users/95656/Desktop/hyp/TZX/")
setwd(path1)
# 计算基因length

gffpath = "Gmax_275_Wm82.a2.v1.gene.gff3"
print(gffpath)
txdb <- makeTxDbFromGFF(gffpath,format="gff")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
length1 = t(as.data.frame(exons_gene_lens))
write.table(length1,"C:/Users/95656/Desktop/hyp/TZX/length.txt", sep = "\t")

# 合并 reads count 数据和length
length1 = read.table("C:/Users/95656/Desktop/hyp/TZX/length.txt", header = T,sep="\t")
countsfile = "C:/Users/95656/Desktop/hyp/TZX/PRJNA238493_BP_raw_count.expression.tsv"
count1 = read.table(countsfile, header = T, sep = "\t")
a = inner_join(count1, length1, by="Locus_Name")
a[is.na(a)] = 0 # 设置 NA值为0
a = a[, -28]  # 删除第28列  Locus_Name.1
write.table(a,"C:/Users/95656/Desktop/hyp/TZX/countsdata.txt", sep = "\t", row.names = FALSE)

# 按照对应关系给countdata columns 改名，改好的文件保存为 re_countsdata.txt
# {'SAMN02644800': 'root', 'SAMN02644811': 'leaf3', 'SAMN02644826': 'seed4', 'SAMN02644825': 'seed3', 'SAMN02644822': 'podseed3', 'SAMN02644821': 'podseed2', 'SAMN02644820': 'podseed1', 'SAMN02644803': 'stem1', 'SAMN02644806': 'leafbud1', 'SAMN02644804': 'stem2', 'SAMN02644813': 'flower2', 'SAMN02644814': 'flower3', 'SAMN02644816': 'flower5', 'SAMN02644815': 'flower4', 'SAMN02644817': 'pod1', 'SAMN02644823': 'seed1', 'SAMN02644818': 'pod2', 'SAMN02644819': 'pod3', 'SAMN02644824': 'seed2', 'SAMN02644809': 'leaf1', 'SAMN02644807': 'leafbud2', 'SAMN02644802': 'cotyledon2', 'SAMN02644810': 'leaf2', 'SAMN02644805': 'shoot meristem', 'SAMN02644812': 'flower1', 'SAMN02644827': 'seed5'}

# 计算 TPM 和 FPKM

setwd("C:/Users/95656/Desktop/hyp/TZX/")
countfile = "C:/Users/95656/Desktop/hyp/TZX/re_countsdata.txt"
mycounts<-read.table(countfile, header = T,sep="\t")
head(mycounts)
rownames(mycounts)<-mycounts[,1]
mycounts<-mycounts[,-1]
head(mycounts)
# TPM
kb <- mycounts$length / 1000
countdata <- mycounts[,1:26]
rpk <- countdata / kb
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)
write.table(tpm,file="my_TPM.txt",sep="\t",quote=F)
# FPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
head(fpkm)
write.table(fpkm,file="my_FPKM.txt",sep="\t",quote=F)


## FPKM转化为TPM
# fpkm_to_tpm = t(t(fpkm)/colSums(fpkm))*10^6
# head(fpkm_to_tpm)
计算 共表达相关性
Text Only
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
path1 = "C:\\Users\\95656\\Desktop\\hyp\\比较"
setwd(path1)
rowdata <- read.table(choose.files(),header=T,sep="\t")
# exp1 = rowdata[,5:]
# "cotyledon1", "leafbud3",
exp1 = rowdata[c("cotyledon2", "flower1", "flower2", "flower3", "flower4", "flower5", "leaf1", "leaf2", "leaf3", "leafbud1", "leafbud2", "pod1", "pod2", "pod3", "podseed1", "podseed2", "podseed3", "root", "seed1", "seed2", "seed3", "seed4", "seed5", "shootmeristem", "stem1", "stem2")]
rownames(exp1)<-rowdata$Name
exp1[is.na(exp1)] = 0 # 设置 NA值为0
exp1 = as.matrix(exp1)
geneid <- rowdata[,'Name']
TCP = geneid[1:5] 
Gma = geneid[1:5]
outfile = choose.files()
for(i in 1:length(TCP))
{
    for(j in 1:length(Gma))
    {
        x <- as.numeric(exp1[i,])
        y <- as.numeric(exp1[j,])
        if(sum(exp1[i,])*sum(exp1[j,]) != 0)#取出表达值都为0的行，如果要严格一点，想要把表达值出现0的行都去掉的话可以把求和函数sum()换成求积函数prod(x)
        {
            gene_cor <- cor.test(x,y,alternative = "two.sided",method = "pearson")
            gene_p <- gene_cor$p.value
            gene_e <- gene_cor$estimate
            if(abs(gene_e) >= 0 && gene_p <= 1)
            {
                res<-paste(TCP[i],Gma[j],gene_e,gene_p,sep="\t")
                write.table(res, outfile,sep="\t",col.names=F,row.names=F,quote=F,append=T)
            }
        }else{
            print("something is error")
        }
    }
}
用 python df 计算相关性
Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import pandas as pd
from scipy import stats
import numpy as np


def get_dfcorr(df1, df2):
    """
    按行计算df1(gene1) 和 df2(gene2) 中数据的相关性, 其中 df1、2的列数可以不相等，method="pearson"
    100W 1m32s
    df 中 列为组织名称，行为基因名称
    """
    df1 = df1.loc[~(df1 == 0).all(axis=1), :]
    df2 = df2.loc[~(df2 == 0).all(axis=1), :]
    df1index = df1.index.tolist()
    df2index = df2.index.tolist()
    list1, list2 = df1.values.tolist(), df2.values.tolist()
    list_r = []
    for n, x in enumerate(list1):
        for m, y in enumerate(list2):
            a = np.array(x).astype(float)
            b = np.array(y).astype(float)
            r, p = stats.pearsonr(a, b)
            listtmp = [df1index[n], df2index[m], r, p]
            list_r.append(listtmp)
    df = pd.DataFrame(list_r, columns=["gene1", "gene2", "correlation", "pvalue"])
    return df
  • 分类于 生物信息学
  • 需要 5 分钟阅读时间

生物信息学常用的Linux命令

1、 前言

东西有点杂乱，还需要慢慢整理。其中的 sed awk perl 等命令我觉得都可以单独分出来讲了

2、 构建blast数据库

Bash
1
2
3
4
makeblastdb -in ../ATaa.fa -dbtype prot -title ATaa-parse_seqids -out ATdb   #可用的蛋白库
makeblastdb -in ../ATaa.fa -dbtype nucl -title ATaa-parse_seqids -out ATdb  #可用的核酸白库
blastp -db /home/hanyapeng/prot/DB/Mtrdb -query At.fa -outfmt 6 -evalue 1e-20 -out Mtr.blast -num_threads 8
blastn -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8

3、 碱基/蛋白序列比对

Text Only
1
2
3
muscle -in AtLj.fa -out AtLj.phy  # -out，输出默认的fasta格式
muscle -in At.fasta -phyiout At.phy  # -phyiout，则输出sequential phylip格式
mafft aa.fa > aa.fas

4、 画树

Bash
1
2
3
nohup raxmlHPC-PTHREADS -s ALO.phy -n ALO -m PROTCATJTT -f a -N 100 -p 12345 -x 1234567 -T 8 &  # RaxML树
nohup ../software/PhyML_3.0_linux64 -i AGO.phy -d aa -b 0  -f e -s BEST -c 4 -m VT -v e -a e -o lr --no_memory_check &  # Phyml 树
../software/FastTree TCP123.phy  > TCP.tree   #fasttree 画树

5、 提取CDS/pep序列

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
nohup hmmsearch -E 1e-5 TCP.hmm  /home/hanyapeng/prot/Gma.prot > Gma.hmm &
less -S Gma.hmm |awk '{if ($1~/e-/) print $9}' > list     #if第一行匹配到e-，则调出第9行#  

nohup hmmsearch --domtblout hmmresult.txt -E 1e-5 TCP.hmm  /home/hanyapeng/prot/Gma.prot  &
less -S hmmresult.txt |awk '{ if ($1~/e-/) print$1"\t"$20"\t"$21}'|less -S > GmaTCPdomain.list    # 提取第1.20.21列为坐标

perl /home/lixiangyong/software/getseq.pl Osa.list /home/hanyapeng/prot/Osa.prot > Osa.fa
perl /home/lixiangyong/software/get_domain.pl  /home/hanyapeng/prot/Osa.prot  OsaTCPdomain.list > OsaTCPdomain.fa  #提取domain

perl /home/lixiangyong/software/getseq.pl Chr13 Gmax275.fa > Chr13.fa  
perl /home/hanyapeng/comm/getpromoter.pl /data/plant_genome/Glyma2.0/Gmax275.fa Q1.txt > 1  # 提取启动子序列 Q1为启动子坐标
cat AT.fa Lj.fa Osa.fa > ATLjOsa.fa

6、 修改环境变量

Text Only
1
2
3
4
5
6
vi .bash_profile     # 或则 .bash_login   .profile   .bashrc 
vi .bash_profile     # 添加环境变量 PATH=$PATH:$HOME/bin
source .bash_profile  #   立刻执行更改

nano .bashrc             #和vi类似的编辑器
source .bashrc

7、 prottest的使用

Text Only
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
nohup java -jar prottest-3.4.jar -i AGO.phy -all-matrices -all-distributions  > a.txt  &  # 此.phy 文件为比对过的氨基酸序列文件

nohup java -jar prottest-3.4.jar -i namefail.phy -o namefail.txt -S 1 -AIC -I -G -IG -F -ncat 8 -all &

nohup java -jar prottest-3.4.jar -i AGOdomain.phy -o AGOdomain.txt -all-matrices -all-distributions &

nohup java -jar jModelTest.jar -d example-data/aP6.fas -f -i -g 4 -s 11 -AIC -a >a.txt  &  此.fas文件为核苷酸序列文件

PhyML-3.1_linux64 -i aP6.phy -d nt -n 1 -b 0 --run_id F81+I -m 000000 -f m -v e -c 1 --no_memory_check -o tlr -s BEST

PhyML-3.1_linux64 -i aP6.phy -q -d aa -m JTT -c 4 -a e

8、 shell 脚本知识

Bash
1
!= 不等于,如:if [ "$a" != "$b" ]

9、 awk的使用

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
less Y2.blast |awk ' {if ($11 ~ /0\.0/) print $0}' > 1      # 每一行中，如果第11列=0.0，则输出全部  “~”是 =  ，$11 调取第11列 ， \ 反义符 
awk '{line[NR]=$0}END{for(i=NR;i>0;i--)print line[i]}' sed.txt  # 把sed.txt 中的内容按行倒着输出
less -S Gma3.fa |awk '{if ($1~/>/) print $head}' > a
awk '!seen[$0]++' <filename>    # awk 去除重复的行
less Step1.2.sh | grep "=" | cut -d "=" -f 1 | awk '{printf $quot" ";$0}' #将一列输出为一行
awk '{ print $(NF-1)}' file     # 出输出倒数第二列

less file |  awk  '{split($1,a,":");print a[1]"-"a[3]"\t"$5}' > our.w.d    # 对awk截取字符串，file中的数据格式为:`ACAN:ENST00000561243:p.D1390E The mutant peptide: APGVEEISGLP is found in wild type sequence`， [参考](https://www.cnblogs.com/sunada2005/p/3493941.html)

awk 'BEGIN{total=0}{total+=$1}END{print total}' file   # 对文件第一列求和

10、 sed的使用

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
sed -i 's/Gma/Glyma/g' Glyma.fa
sed -i '/^$/d' 2_domain.txt  # 删除文件空行
sed -i '/#/d' AtTCP.txt      # 删除有\#的行
sed -i s/^M// AAt.list.txt   # 此处^M必需手动输入！！
sed -i /^$/d file     #删除空行   d 指删除含有空行符的这一行
sed -i s/tr.*//g LjTCP.fa  #  去除无用的注释    :%s/|.*.\[//g  #去除从|开始，到[结束内的所有字符，包括[#

sed 's/^/aaa/g' filename.txt > outputfile.txt  # 在文件中的每一行之前加上一个字符串，比如：aaa    行尾将^是 $ 
sed -i/\>Chr17/d Chr17.fa   # 删除Chr17.fa 中的Chr字符
sed -i "s/zhangsan/lisi/g" `grep zhangsan -rl modules` #批量替换，讲文件夹modules中所有文件执行sed -i       "s/zhangsan/lisi/g" 该符号`不是点，是Esc下面的符号 
sed -i '1iNew Top Line' aa.txt    #  sed 在句首添加一行 New Top Line

11、 linux 命令

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
date   #  查看系统日期
cal 1 2016   # 查看2016年1月的日历
df     # 查看磁盘驱动器的可用空间  
free   # 显示可用内存
.     # 工作目录     
..   # 工作目录的父目录   如：  cd ..   返回上一级    cd  ./software/  
file Gma1.fa   # 显示  Gma1.fa: ASCII text  该文件为文本文档ASCII
alias foo='cd /usr; ls; cd -'   # 自己创建命令   unalias foo 删除
gedit Gma1.fa   # 图形编辑器编辑Gma1.fa 

rm -- -foo   # 要删除第一个字符为‘-’的文件 (例如‘-foo’)
rm ./-foo   # 要删除第一个字符为‘-’的文件 (例如‘-foo’)

chmod 755 lixiangyong/  # 更改目录权限为可读

wc -l filename # 就是查看文件里有多少行
wc -w filename # 看文件里有多少个word。
wc -L filename # 文件里最长的那一行是多少个字
wc -c  filename   # 统计字节数

cat Y.list | sort -u > YRedundancy.list    # 去除文件Y.list中的重复的行命名为YR.list
less -S AGO.fa |grep ">" |wc -l  # 检索>有多少个  同  grep -c ">" AGO.fa
for tar in *.tar.gz;  do tar xvf $tar; done    #批量解压文件#
split -l 10 splitfile.txt D  # 将文件按每10行分割，新文件命名为 Daa Dab#
rename linux unix linux*    # 文件批量  
find . -name '*.txt' | xargs perl -pi -e 's/aaa/bbb/g'  # 遍历整个文件夹，把所有txt文件中的aaa替换成bbb
mkdir $(date +%F)  #创建一个以日期为名的文件夹
touch {1..9}.txt   # 批量创建文件

ls -F | grep '/$'         #  查看当前目录下的文件夹 
for i in `ls -F | grep '/$'`;do sudo du -sh $i;done   #查看当前目录下，每个子目录及其子目录的总大小

rename "s/AA/aa/" *
rename .fas .fa *.fas    # CentOS 下可用

iostat -d -k 1 10可以看到每个盘的读写速度

find /home/lxy -type f -size +100M    #列出文件夹 /home/lxy 下大于100M的文件
find . -type l -exec rm  {} \;  #列出当前路径 .  下的link类型的文件(type
 l) 并删除

12、 单行 perl

Text Only
1
perl -e "print lc($a)"    # Perl 大小写转换

13、 R 命令行

Text Only
1
Rscript /home/muscle/test.r  # linux 上运行r脚本

14、 kaks计算

Text Only
1
2
perl pal2nal.perl aligned.protein.fas original.cds.fas -output fasta > aligned.cds.fas
perl ../software/pal2nal.pl 2.fas 2.cds -output paml > 2.nuc    #kaks计算

15、 gff和gtf之间的格式转化

Text Only
1
2
gffread my.gff3 -T -o my.gtf # gff2gtf
gffread merged.gtf -o- > merged.gff3  # gtf2gff

16、 vim 的使用

Text Only
1
2
3
4
5
6
7
8
9
:12,23s/aa/bb/g     # 将从12 行到23 行中出现的所有包含aa 的字符串中的aa 替换为bb
:s/\<aa\>/bb/g      # 将光标所在行出现的所有aa 替换为bb, 仅替换aa 这个单词
:12,23s/^/#/        # 将从12 行到23 行的行首加入# 字符

:10,15s#^#\#\##g          #  vim的 V模式下，下注释x  10 -15  行，##注释

Esc shift+ : wq   enter   =  Esc shift+ : x  enter       #保存并退出
Esc shift+ : q！  enter                                  #不保存强制退出 
vi foo.txt ls-output.txt    # 同时vi编辑多个文件  ：n   ：N   ：buffer 2 来回切换     也可以先打开一个，在用：e ls-output.txt

17、 数据存放

/data/public/hanyapeng/Gm/Gma.collinear.groups # 大豆中所有旁系同源基因
/data/public/wanglei/inparanoid4.0/Soybeanparalog2.txt # 大豆中所有旁系同源基因

推荐连接