Pytorch是目前最火的深度学习框架之一，另一个是TensorFlow。不过我之前一直用到是CPU版本，几个月前买了一台3070Ti的笔记本（是的，我在40系显卡出来的时候，买了30系，这确实一言难尽），同时我也有一台M1芯片Macbook Pro，目前也支持了pytorch的GPU加速，所以我就想着，在这两个电脑上装个Pytorch，浅度学习深度学习。

Apple silicon

首先是M1芯片，这个就特别简单了。先装一个conda，只不过是内置mamba包管理器，添加conda-forge频道，arm64版本。

1
2
3
4

# 下载
wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-MacOSX-arm64.sh
# 安装
bash Mambaforge-MacOSX-arm64.sh

然后我们用mamba创建一个环境，用的是开发版的pytorch，所以频道指定pytorch-nightly

1
2
3

mamba create -n pytorch \
   jupyterlab jupyterhub pytorch torchvision torchaudio 
   -c pytorch-nightly

最后，用conda activate pytorch,然后测试是否正确识别到GPU

1
2
3
4
5

import torch
torch.has_mps
# True
# 配置device
device = torch.device("mps")

参考资料: https://developer.apple.com/metal/pytorch/

Windows NVIDIA

首先，需要确保你的电脑安装的是NVIDIA的显卡，以及有了相应的CUDA驱动。

CUDA的显卡架构要求: https://docs.nvidia.com/deeplearning/cudnn/support-matrix/index.html

新一代的电脑上基本都自带CUDA驱动。可以通过打开NVIDIA控制面板的系统信息

在组件中查看你已经安装的CUDA驱动，例如我的是11.7.89 。

也可以通过命令行的方式查看，

接下来，我们来安装pytorch。同样也是推荐conda的方法，我们先从清华镜像源中下载Miniconda。

选择Windows的安装包

安装完之后，我们就可以通过Anaconda Prompt进入命令行，根据pytorch网站上的推荐进行安装。

但是有一点不同，为了避免环境冲突，最好是单独创建一个环境，所以代码如下

1

conda create -n pytorch pytorch torchvision torchaudio pytorch-cuda=11.7 -c pytorch -c nvidia

接着用 conda activate pytorch启动环境，然后在python环境中测试

1
2
3

import torch
torch.has_cuda
# True

几个常见的问题（至少我在写文章想到的问题）：

Q: 使用conda和pip安装的区别是什么？

A: conda是pytorch官方推荐的安装方式，因为conda会一并帮你装好pytorch运行所需要的CUDA驱动和相关工具集。这意味着为conda所占用的空间会更多一些。

Q: 既然conda那么强，是不是让conda帮我在非常老的硬件上安装最新的pytorch吗？

A: 我觉得这个跟装游戏类似，你虽然能装上游戏，但是不满足游戏的最低配置需求，照样跑不动。放在conda上，conda虽然能给你配置好cuda驱动，但是cuda驱动本身对系统的显卡驱动有要求。如果你的显卡驱动不满足，那么就算装好了，也会报错。RuntimeError: The NVIDIA driver on your system is too old

Q: 电脑上必须要安装CUDA驱动和安装CUDA toolkit吗？

A: 其实我个人不是很确定如何回答，如下是我目前的一些见解。如果你用的conda，那么他会帮你解决一些依赖问题。如果你是用pip，那么就需要你自己动手配置。其中，CUDA驱动是必须要安装的，因为CUDA驱动负责将GPU硬件与计算机操作系统相连接，不装驱动，操作系统就不识别CUDA核心，相当你没装NVIDIA显卡。而CUDA toolkit是方便我们调用CUDA核心的各种开发工具集合，你装CUDA toolkit的同时会配套安装CUDA驱动。除非你要做底层开发，或者你需要从源码编译一个pytorch，否则我们大可不装CUDA toolkit。

Q: 如果我电脑上的CUDA驱动版本比较旧怎么办？或者说我CUDA的驱动版本是11.7，但是我安装了cuda=11.8的pytorch，或者版本不一样的pytorch会怎么样？

A: 我们在安装cuda=11.7的pytorch，本质上安装的是在CUDA Toolkit版本为11.7环境下编译好的pytorch版本。当cuda版本之间的差异不是特别的大，或者说不是破坏性的升级的情况下，那么理论上也是能运行的。

手写数据性能测试

下面用的是GPT3.5给我提供一段手写字识别（MNIST)案例代码，用来测试不同的平台的速度。

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

import torch
import torchvision
import torchvision.transforms as transforms

# 转换数据格式并且加载数据
transform = transforms.Compose(
    [transforms.ToTensor(),
     transforms.Normalize((0.5,), (0.5,))])

trainset = torchvision.datasets.MNIST(root='./data', train=True,
                                        download=True, transform=transform)
trainloader = torch.utils.data.DataLoader(trainset, batch_size=64,
                                          shuffle=True, num_workers=2)

testset = torchvision.datasets.MNIST(root='./data', train=False,
                                       download=False, transform=transform)
testloader = torch.utils.data.DataLoader(testset, batch_size=64,
                                         shuffle=False, num_workers=2)

# 定义网络模型
class Net(torch.nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv1 = torch.nn.Conv2d(1, 6, 5)
        self.pool = torch.nn.MaxPool2d(2, 2)
        self.conv2 = torch.nn.Conv2d(6, 16, 5)
        self.fc1 = torch.nn.Linear(16 * 4 * 4, 120)
        self.fc2 = torch.nn.Linear(120, 84)
        self.fc3 = torch.nn.Linear(84, 10)

    def forward(self, x):
        x = self.pool(torch.nn.functional.relu(self.conv1(x)))
        x = self.pool(torch.nn.functional.relu(self.conv2(x)))
        x = x.view(-1, 16 * 4 * 4)
        x = torch.nn.functional.relu(self.fc1(x))
        x = torch.nn.functional.relu(self.fc2(x))
        x = self.fc3(x)
        return x

net = Net()

# 这里的代码比较随意，就是用哪个平台运行哪个
# CPU
device = torch.device("cpu")
# CUDA
device = torch.device("cuda:0")
# MPS
device = torch.device("mps")

net.to(device)

# 定义损失函数和优化器
criterion = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.SGD(net.parameters(), lr=0.001, momentum=0.9)

# 训练网络

import time

start_time = time.time()  # 记录开始时间

for epoch in range(10):  # 进行10次迭代训练
    running_loss = 0.0
    for i, data in enumerate(trainloader, 0):
        inputs, labels = data[0].to(device), data[1].to(device)
        optimizer.zero_grad()
        outputs = net(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()

        running_loss += loss.item()
        if i % 100 == 99:
            print('[%d, %5d] loss: %.3f' %
                  (epoch + 1, i + 1, running_loss / 100))
            running_loss = 0.0

end_time = time.time()  # 记录结束时间
training_time = end_time - start_time  # 计算训练时间

print('Training took %.2f seconds.' % training_time)

print('Finished Training')

# 测试网络
correct = 0
total = 0
with torch.no_grad():
    for data in testloader:
        images, labels = data[0].to(device), data[1].to(device)
        outputs = net(images)
        _, predicted = torch.max(outputs.data, 1)
        total += labels.size(0)
        correct += (predicted == labels).sum().item()

print('Accuracy of the network on the 10000 test images: %d %%' % (
      ))

最后统计下来的结果为

Windows平台

• 3070Ti Training took 45.02 seconds.
• i9 12900H Training took CPU 75.65

Mac平台

• M1 Max Training took 50.79 seconds.
• M1 MAX CPU Training took 109.61 seconds.

总体上来说，GPU加速很明显，无论是mac还是Windows。其次是GPU加速效果对比，M1 max芯片比3070Ti差个10%？。

不过目前测试都是小数据集，等我学一段时间，试试大数据集的效果。

先提一个问题，请问下面最后的两个逻辑语句的判断，结果是什么？

1
2
3
4
5
6
7
8
9
10

x <- c(1234567.8, 12.12345)
x11 <- format(x, scientific = FALSE)
x12 <- format(x, scientific = TRUE)
x21 <- as.numeric(x11)
x22 <- as.numeric(x12)

x21 == x
x22 == x

大家可以稍作思考。

之所以给出开头的问题，还是因为之前讨论的apply函数。当apply应用到data.frame上时，中间会调用泛型函数as.matrix将数据进行转换，如果发现数据中存在非数值型，会调用format对数值型数据进行格式化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

....
  else if (non.numeric) {
    for (j in pseq) {
      # 如果是字符串，不转换
      if (is.character(X[[j]])) 
        next
      # 如果是逻辑值，调用as.character
      else if (is.logical(xj <- X[[j]])) 
        xj <- as.character(xj)
      else {
        # 如果有NA，先确定位置
        miss <- is.na(xj)
        # 如果是因子，调用as.vector，否则调用format
        xj <- if (length(levels(xj))) 
          as.vector(xj)
        else format(xj)
        # 最后把NA放回
        is.na(xj) <- miss
      }
      X[[j]] <- xj
    }
  }
....

这个format的描述只有一句，‘Format an R object for pretty printing.’即，格式化一个 R 对象用于漂亮的打印，

问题出在这个“漂亮”上，R语言有自己的理解。

比如说同样的 1234567.8在不同的环境下，会出现了不同的情况，

1
2
3
4
5

format(c(1234567.8, 12.12345))
# "1.234568e+06" "1.212345e+01"
format(c(1234567.8, 123.12345))
# "1234567.8000" "    123.1235"

一个采取了科学计数法，并保留了7位有效数字，显示为1.234568e+06，一个则保持原状1234567.8000。

为什么会出现这个情况呢？这和scientific参数有关，接受一个逻辑值或者一个整数。如果提供的逻辑值，行为非常简单，就是要么用科学计数法，要么不用科学计数法。

而问题出在它默认值不是逻辑值，而是NA，此时他使用 getOption("scipen") 的结果作为输入，用于决定是以固定，还是以指数符号（科学计数法）打印数字。正值偏向于固定符号，负值则偏向于科学计数法。

具体这个参数会触发R哪一行底代码呢？通过我刨根问底，最终确定源代码中src/mian/format.c 的formatReal函数，一个用于对double类型数据进行格式化的函数。

代码有一段注释，含义是，在满足有效位数的前提下，固定表示使用的空间不超过指数表示时，优先使用固定显示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

    /* F Format: use "F" format WHENEVER we use not more space than 'E'
     *    and still satisfy 'R_print.digits' {but as if nsmall==0 !}
     *
     * E Format has the form   [S]X[.XXX]E+XX[X]
     *
     * This is indicated by setting *e to non-zero (usually 1)
     * If the additional exponent digit is required *e is set to 2
     */
     ...
  *d = mxns - 1;
  *w = neg + (*d > 0) + *d + 4 + *e; /* width for E format */
  if (wF <= *w + R_print.scipen) { /* Fixpoint if it needs less space */
      *e = 0;
      if (nsmall > rgt) {
    rgt = nsmall;
    wF = mxsl + rgt + (rgt != 0);
      }
      *d = rgt;
      *w = wF;
  } /* else : "E" Exponential format -- all done above */
...       

代码中w是指数表示法的宽度, wF是固定表示的宽度，最后w会加上 scipen 得到最终的长度，默认是0，也就是公平决斗。

回到数据c(1234567.8, 12.12345)中，我们来看它在固定表示和指数表示的结果。

1
2
3
4
5
6
7

> format(c(1234567.8, 12.12345), scientific = TRUE)
[1] "1.234568e+06"
[2] "1.212345e+01"
> format(c(1234567.8, 12.12345), scientific = FALSE)
[1] "1234567.80000"
[2] "     12.12345"

在保证最小数字有7位有效数字的前提下（由digits参数控制），采用指数显示法时结果是12个字符，采取固定显示结果是13个字符。当scipen=0时，指数显示比较窄，因此选择指数显示。

而数据c(1234567.8, 123.12345)，在保证最小数字有7位有效数字时，结果都是12个字符，势均力敌的情况下，采取固定显示。

1
2
3
4
5
6

> format(c(1234567.8, 123.12345), scientific = TRUE)
[1] "1.234568e+06"
[2] "1.231235e+02"
> format(c(1234567.8, 123.12345),scientific = FALSE)
[1] "1234567.8000"
[2] "    123.1235"

有效数字的定义是，从第一位非0数字开始到最后一位，一共有多少位数字。例如100是3为有效数字，0.001是1位有效数字，0.00100是3位有效数字。

最后帮助大家梳理下整体逻辑

1. format函数作用是让打印结果更加美观
2. 在打印数字时，在保证最小数以一定有效数值展示的前提下，对比固定表示和指数表示两者的字符数目，选择占用字符最少的方法
3. 最终数值的展示效果会受到digits和scientific的影响
4. digits和scientific默认值来自于getOption(“digits”)和getOption(“scipen”)
5. 如果想让结果都以固定显示，可以设置scipen为一个非常大的数。

我使用GPT也有一段时间了（主要是3.5，还在等4），简单说说我的使用感觉。

• 你需要把GPT的回答都当做一次闭卷考试，因此它不能当场查资料。它为了得分，需要把每道题目都写出来，无论它有没有把握。我们都应该记得老师说过，要把试卷写满吧。
• 因此，有一些题目他准备的非常好（指训练集里面有非常多的语料），有一些题目他不是特别的熟练（语料适中），而有些他压根不会（没有语料资源）。
• 你不能让他做一些记忆性的题目，因为他背的东西可能太多了，以至于很多内容都背混了。
• 有些时候，他非常的倔强，觉得自己就是对的。

事情起因是这样子，我有个genbank文件，里面记录的是线粒体的注释的坐标信息。由于这个注释是来源同源注释，可能存在错误，因此还需要后续的手动修改。我没去找到专门的工具，只是根据分析结果，手动改了里面一些基因的坐标，就导致里面的translation这一栏跟坐标对不上，需要更新一下。虽然我有能力写这个代码去更新结果，但是我希望给GPT(模型turbo 3.5)一个机会，让他处理下。

首先，我问了它关于genbank的格式定义，确保它后续能够进行处理。他的回答非常正确，符合我的认知（下面是对话截图）

接下来我就开始问问题了，我让他编写脚本，让我调整坐标，它写了一段python 的代码，通过调用BioPython模块来实现我刚才说的功能。

但是这代码有问题，因为通用性比较差，所以我让他改写代码，有一个遍历过程，先得到新的CDS，然后和原来的translation进行比较，如果存在不同则替换

它的工作非常棒，下面是它给的结果。

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

from Bio import SeqIO

# 读取GenBank文件
record = SeqIO.read("your_genbank_file.gb", "genbank")

# 遍历所有CDS区域
for feature in record.features:
    if feature.type == "CDS":
        # 获取原始位置信息
        old_start = feature.location.start
        old_end = feature.location.end
        
        # 修改位置信息为新的位置
        new_start = 100  # 新的起始位置
        new_end = 500  # 新的终止位置
        feature.location = feature.location._shifted(new_start - old_start)
        feature.location = feature.location._shifted(new_end - old_end)
        
        # 更新翻译信息
        old_translation = feature.qualifiers.get("translation", [""])[0]
        new_sequence = feature.extract(record.seq)  # 提取新序列
        new_translation = new_sequence.translate(to_stop=True)  # 计算新的翻译信息
        new_translation_str = str(new_translation)
        if new_translation_str != old_translation:
            feature.qualifiers["translation"] = [new_translation_str]

# 写回GenBank文件
with open("modified_genbank_file.gb", "w") as output_handle:
    SeqIO.write(record, output_handle, "genbank")

不过它总是念念不忘他的new_start 和 new_ned，这显然不是我需要的，所以，我让他把这部分去掉。下面截图是我的提问和它的回答。

这次代码简洁多了。不过他还不是一个合格的代码，因为调用只能是固定的文件名，我希望它能进一步提高通用性，所以我的提问是“好的，进一步提供他的通用性，使得我们可以在命令行通过 python remake.py input.gb output.gb来实现文件的重新生成”

这次他的回答更加的棒，它用了argparse模块，并将原来的代码函数化，下面是他的回答。

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

import argparse
from Bio import SeqIO


def update_cds_translation(record, new_positions):
    """
    更新GenBank记录中所有CDS区域的翻译信息并返回修改后的记录。
    :param record: GenBank记录对象
    :param new_positions: 包含每个CDS区域新起始和终止位置的dict对象
    :return: 修改后的GenBank记录对象
    """
    for feature in record.features:
        if feature.type == "CDS":
            # 获取CDS区域原始序列和翻译信息
            old_sequence = feature.extract(record.seq)
            old_translation = feature.qualifiers.get("translation", [""])[0]

            # 获取CDS区域新位置
            start, end = new_positions.get(feature.qualifiers['gene'][0], (feature.location.start, feature.location.end))

            # 调整CDS区域位置并计算新翻译信息
            feature.location = feature.location._shifted(start - feature.location.start)
            feature.location = feature.location._shifted(end - feature.location.end)
            new_sequence = feature.extract(record.seq)
            new_translation = new_sequence.translate(to_stop=True)
            new_translation_str = str(new_translation)

            # 如果翻译信息不同，则更新qualifiers字段
            if new_translation_str != old_translation:
                feature.qualifiers["translation"] = [new_translation_str]

    return record


if __name__ == '__main__':
    # 解析命令行参数
    parser = argparse.ArgumentParser(description="Update CDS translations in a GenBank file.")
    parser.add_argument("input", help="Input GenBank file name")
    parser.add_argument("output", help="Output GenBank file name")
    parser.add_argument("-p", "--positions", help="New start and end positions for each gene, e.g. 'gene1:100-500 gene2:200-600'", default="")
    args = parser.parse_args()

    # 读取GenBank文件和新位置信息
    record = SeqIO.read(args.input, "genbank")
    new_positions = {gene_pos.split(":")[0]: tuple(map(int, gene_pos.split(":")[1].split("-"))) for gene_pos in args.positions.strip().split()}

    # 更新CDS翻译信息并写回GenBank文件
    updated_record = update_cds_translation(record, new_positions)
    with open(args.output, "w") as output_handle:
        SeqIO.write(updated_record, output_handle, "genbank")

更惊喜的是，在它的解释中，它还扩展了功能，它说道

此代码使用argparse模块解析命令行参数，并根据提供的输入文件名和输出文件名读取和写入GenBank记录。它还允许通过命令行参数-p/–positions指定每个基因的新起始和终止位置，方式为’gene1:100-500 gene2:200-600’。如果不提供此参数，则使用每个CDS区域的原始位置。最后，它调用update_cds_translation函数来更新所有CDS区域的翻译信息，然后将修改后的GenBank记录写回到输出文件中。

我觉得差不多了，并复制了代码，开始运行，然后我报错了。错误信息如下

1
2
3
4
5
6
7

我运行了代码，出现了如下报错，我应该如何处理
Traceback (most recent call last):
  File "/Users/xuzhougeng/Downloads/remake.py", line 48, in <module>
    updated_record = update_cds_translation(record, new_positions)
  File "/Users/xuzhougeng/Downloads/remake.py", line 22, in update_cds_translation
    feature.location = feature.location._shifted(start - feature.location.start)
AttributeError: 'FeatureLocation' object has no attribute '_shifted'

我直接复制了这段代码，问gpt怎么办，它推断是biopython的版本问题，还给我写了一个更加通用的代码，新的代码如下

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

import argparse
from Bio import SeqIO, SeqFeature


def update_cds_translation(record, new_positions):
    """
    更新GenBank记录中所有CDS区域的翻译信息并返回修改后的记录。
    :param record: GenBank记录对象
    :param new_positions: 包含每个CDS区域新起始和终止位置的dict对象
    :return: 修改后的GenBank记录对象
    """
    for feature in record.features:
        if feature.type == "CDS":
            # 获取CDS区域原始序列和翻译信息
            old_sequence = feature.extract(record.seq)
            old_translation = feature.qualifiers.get("translation", [""])[0]

            # 获取CDS区域新位置
            start, end = new_positions.get(feature.qualifiers['gene'][0], (feature.location.start, feature.location.end))

            # 调整CDS区域位置并计算新翻译信息
            offset_start = start - feature.location.start
            offset_end = end - feature.location.end
            new_location = SeqFeature.FeatureLocation(
                feature.location.start + offset_start,
                feature.location.end + offset_end,
                feature.location.strand
            )
            new_sequence = feature.extract(record.seq)
            new_translation = new_sequence.translate(to_stop=True)
            new_translation_str = str(new_translation)

            # 如果翻译信息不同，则更新qualifiers字段
            if new_translation_str != old_translation:
                feature.qualifiers["translation"] = [new_translation_str]

            # 更新位置信息
            feature.location = new_location

    return record


if __name__ == '__main__':
    # 解析命令行参数
    parser = argparse.ArgumentParser(description="Update CDS translations in a GenBank file.")
    parser.add_argument("input", help="Input GenBank file name")
    parser.add_argument("output", help="Output GenBank file name")
    parser.add_argument("-p", "--positions", help="New start and end positions for each gene, e.g. 'gene1:100-500 gene2:200-600'", default="")
    args = parser.parse_args()

    # 读取GenBank文件和新位置信息
    record = SeqIO.read(args.input, "genbank")
    new_positions = {gene_pos.split(":")[0]: tuple(map(int, gene_pos.split(":")[1].split("-"))) for gene_pos in args.positions.strip().split()}

    # 更新CDS翻译信息并写回GenBank文件
    updated_record = update_cds_translation(record, new_positions)
    with open(args.output, "w") as output_handle:
        SeqIO.write(updated_record, output_handle, "genbank")

这个代码对不对呢？至少我跑完是没有报错的，我进一步的检查了它的分析结果。发现很多CDS的长度明显短了，我觉得这可能是我没有说清楚，我给的是线粒体基因组。我讲这个结果反馈给了gpt，它认可了我的推测，下面是他的回答。

后来我们又反反复复商量了好几轮，对代码反复打磨，通常都是它给出了建议，我询问他是不是有另一种解决思路，比如说有些genbank文件里的cds长度有问题，不一定是3的倍数，他建议我对处理来的文件，我觉得可以考虑捕获报错，他也执行的特别好。

下面的代码目前我最为满意的版本了

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
68

import argparse
from Bio import SeqIO, Data


def update_cds_translation(record, new_positions):
    """
    更新GenBank记录中所有CDS区域的翻译信息并返回修改后的记录。
    :param record: GenBank记录对象
    :param new_positions: 包含每个CDS区域新起始和终止位置的dict对象
    :return: 修改后的GenBank记录对象
    """
    # 遍历每个特征
    for feature in record.features:
        # 如果是CDS区域
        if feature.type == 'CDS':
            # 获取特征的名称
            gene = feature.qualifiers.get('gene', [''])[0]
            # 获取CDS区域的起始和终止位置
            start, end = feature.location.start, feature.location.end
            # 如果有新位置信息，则使用新位置信息
            if gene in new_positions:
                start, end = new_positions[gene]
            # 获取CDS区域的序列
            cds_seq = record.seq[start:end]
            # 判断序列的方向
            if feature.strand == -1:
                cds_seq = cds_seq.reverse_complement()

            # 获取codon_start和transl_table信息
            codon_start = int(feature.qualifiers.get('codon_start', [1])[0])
            transl_table = int(feature.qualifiers.get('transl_table', [1])[0])
            # 翻译CDS区域
            table = Data.CodonTable.unambiguous_dna_by_id[transl_table]
            try:
                cds_protein = cds_seq[codon_start - 1:].translate(table=table, cds=True)
            except Data.CodonTable.TranslationError:
                cds_protein = ''
                print(f"Warning: Unable to translate CDS region {gene} in record {record.id}")
            # 更新CDS区域的翻译信息
            feature.qualifiers['translation'] = [str(cds_protein)]
    return record


if __name__ == '__main__':
    # 解析命令行参数
    parser = argparse.ArgumentParser(description="Update CDS translations in a GenBank file.")
    parser.add_argument("input", help="Input GenBank file name", type=argparse.FileType('r'))
    parser.add_argument("output", help="Output GenBank file name", type=argparse.FileType('w'))
    parser.add_argument("-p", "--positions", help="New start and end positions for each gene, e.g. 'gene1:100-500 gene2:200-600'", default="")
    args = parser.parse_args()

    # 读取GenBank文件和新位置信息
    records = (record for record in SeqIO.parse(args.input, "genbank"))
    new_positions = {gene: tuple(map(int, pos.split('-'))) for gene, _, pos in map(str.partition, args.positions.split())}
    
    # 更新CDS翻译信息并写回GenBank文件
    untranslated_cds_count = 0
    for record in records:
        record = update_cds_translation(record, new_positions)
        # 统计无法翻译的CDS区域个数
        for feature in record.features:
            if feature.type == 'CDS' and not feature.qualifiers.get('translation', [''])[0]:
                untranslated_cds_count += 1
        SeqIO.write(record, args.output, "genbank")
    
    # 提示用户有多少个CDS区域无法翻译为蛋白质序列
    if untranslated_cds_count > 0:
        print(f"Warning: Failed to translate {untranslated_cds_count} CDS regions.")

我运行了这个代码，发现结果非常好。

1
2
3

$ python remake.py mitoscaf.fa.gbf new.gff
Warning: Unable to translate CDS region COX2 in record XZG
Warning: Failed to translate 1 CDS regions.

其实早在几年前，openAI推出GPT3的时候，我就知道了，我也提交了申请，但是一直没被通过，也就被我忘了。后来chatGPT横空出世，我开始也聊了几句，后面也就放在一边了。直到最近，它又火了，并且还出现了一个非常优秀的github项目，chatgpt-web，使得我能搭建了一个自己的专属GPT，之后我还搞定了付费API。

虽然他也会出错（未来的4.0 API或许会好很多），但这样子说的仿佛人不出错一样。我现在就把它当做一个合作搭档，我有一个想法，让我帮我快速实现一下，然后我根据结果进行修改，并反馈给他，他也能迅速的给我一个反馈，这不挺好吗？

利用后端的SpringBoot框架和前端的React框架实现一个文件上传的小项目，对目前的学习做一个总结。

环境搭建

我们先分别搭建前后端的项目环境，然后进行项目开发。

前端环境搭建

前端开发需要一些基本的HTML, CSS和Javscript的背景知识。

我们直接使用成前端比较热门的React框架进行网页搭建，还有一个热门框架是Vue。

先安装NODE环境，在https://nodejs.org/en/下载长期支持版，安装之后，在命令行输入node -v 确定版本。

编程环境选择VSCode，同时安装了几个插件：

• IntelliJ IDEA Keybindings
• Prettier
• Prettier ESlint

之后，我们使用umi作为项目的脚手架，搭建项目

1
2
3
4
5
6

mkdir file_upload_frontend && cd file_upload_frontend
pnpm dlx create-umi@latest
# 选择simple App + pnpm + taobao
# 安装ant-design
pnpm i antd @ant-design/icons

安装结束后，启动项目，访问Network里的地址，就确定项目是跑上了。

1
2
3
4
5
6
7
8
9
10
11

pnpm dev
....
        ╔════════════════════════════════════════════════════╗
        ║ App listening at:                                  ║
        ║  >   Local: http://localhost:8000                  ║
ready - ║  > Network: http://192.168.50.97:8000              ║
        ║                                                    ║
        ║ Now you can open browser with the above addresses↑ ║
        ╚════════════════════════════════════════════════════╝
...

网页显示如下

后端环境搭建

后端使用Java SpringBoot实现RESTful的接口开发。(这里的背景知识要求比较多，Java的面向对象编程，反射，注解，Spring框架，SpringMVC框架）

首先，到https://www.jetbrains.com/idea/download/里下载一个社区版的IntelliJ IDEA。相比于商业版，社区办缺少了Spring初始化、SQL数据库连接、HTTP客户端等功能，但是也是能用。

接着，在https://start.spring.io中创建一个SpringBoot项目，配置如下

得到一个fileupload.zip文件，解压缩后就得到了我们项目的基本结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

fileupload
├── HELP.md
├── mvnw
├── mvnw.cmd
├── pom.xml
└── src
    ├── main
    │   ├── java
    │   │   └── top
    │   │       └── xuzhougeng
    │   │           └── fileupload
    │   │               └── FileuploadApplication.java
    │   └── resources
    │       ├── application.properties
    │       ├── static
    │       └── templates
    └── test
        └── java
            └── top
                └── xuzhougeng
                    └── fileupload
                        └── FileuploadApplicationTests.java

利用IDEA打开这个项目，把resources下的application.properties改成application.yaml，编辑内容如下

1
2

server:
  port: 8070

然后点击IDEA的运行键，如果能正常访问localhost:8070就意味着后端环境配置成功

代码编写

前端基础代码

前端的工作非常简单，只需要去Ant Design的官方找一个合适的模版。

打开https://ant-design.gitee.io/components/overview-cn/，找到Upload上传组件块，选择一个顺眼的，复制它的代码

之后在项目中建立一个src/pages/upload.tsx文件，把代码粘贴进去

修改 .umirc.tsx 中的路由配置

1
2
3
4
5
6
7
8
9
10
11

import { defineConfig } from "umi";

export default defineConfig({
  routes: [
    { path: "/", component: "index" },
    { path: "/docs", component: "docs" },
    { path: "/upload", component: "upload"}
  ],
  npmClient: 'pnpm',
});

打开http://localhost:8000/upload 如果看到里面内容跟复制的内容一样，就算成功。

接下来，编辑upload.tsx文件，做一些基本的修改

1
2

multiple: false, //只允许传一个文件
action: 'http://localhost:8070/files', //设置上传的站点

multiple设置为false，用于限制单次上传的文件数，action设置上传的目标地址。

这样子，前端就算完工了。

后端代码

后端文件上传代码参考https://www.bezkoder.com/spring-boot-file-upload/

首先，我们写服务层的代码，实现两个功能

1. 初始化文件上传后存放的目录
2. 保存文件

在top.xuzhougeng.fileupload建立一个软件包, service。

然后在service中创建一个接口FileStorageService，里面就定义两个接口函数

1
2
3
4
5
6
7
8
9
10

// top/xuzhougeng/fileupload/service/FileStorageService.java
package top.xuzhougeng.fileupload.service;
import org.springframework.web.multipart.MultipartFile;

public interface FileStorageService {

    public void init();
    public void save(MultipartFile file);
}

通过在FileStorageService上用快捷键option + enter，实现接口。

实现类存放在service的impl目录下

代码如下

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

// top/xuzhougeng/fileupload/service/FileStorageServiceImpl.java
package top.xuzhougeng.fileupload.service.impl;

import org.springframework.stereotype.Service;
import org.springframework.web.multipart.MultipartFile;
import top.xuzhougeng.fileupload.service.FileStorageService;

import java.io.IOException;
import java.nio.file.FileAlreadyExistsException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.nio.file.Paths;


@Service
public class FileStorageServiceImpl implements FileStorageService {

    private final Path root = Paths.get("uploads");
    @Override
    public void init() {
        try {
            Files.createDirectories(root);
        } catch (IOException e){
            throw new RuntimeException("无法初始化upload的上传目录");
        }

    }

    @Override
    public void save(MultipartFile file) {
        try {
            Files.copy(file.getInputStream(),
                    this.root.resolve(file.getOriginalFilename()));
        } catch (Exception e){
            if (e instanceof FileAlreadyExistsException){
                throw new RuntimeException("文件已存在");
            }

            throw new RuntimeException(e.getMessage());
        }
    }
}

初始化init调用Files模块创建文件夹。save函数的参数类型为MultipartFile，是SpringMVC用来简化文件上传操作的工作类。通过调用Files模块，将上传的文件流复制到指定的目录下。

服务层搞定之后，就是控制层，创建一个软件包controller，并在目录下建立一个FileController类。

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

package top.xuzhougeng.fileupload.controller;


import org.springframework.beans.factory.annotation.Autowired;
import org.springframework.stereotype.Controller;
import org.springframework.web.bind.annotation.*;
import org.springframework.web.multipart.MultipartFile;
import top.xuzhougeng.fileupload.service.FileStorageService;


@RestController
@CrossOrigin(value = {"http://192.168.50.97:8000"})
public class FileController {

    private FileStorageService fileStorageService;
    @Autowired
    public void setFileStorageService(FileStorageService fileStorageService) {
        this.fileStorageService = fileStorageService;
    }

    @PostMapping("/files")
    public String fileSave(@RequestParam("file") MultipartFile file){
        System.out.println("loading file");
        try {
            fileStorageService.save(file);
            return "Upload successfully";
        } catch (Exception e){
            return "Upload failed";
        }

    }
}

它的功能非常简单，就是相应POST请求，然后调用服务层的fileStorageService的save方法进行文件保存。

最后还需要修改启动函数，使其实现CommandLineRunner，覆写他的 run方法，使其在运行前调用FileStorageService的init方法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

@SpringBootApplication
public class FileuploadApplication implements CommandLineRunner {

  @Resource
  private FileStorageService fileStorageService;

  @Override
  public void run(String... args) throws Exception {
    fileStorageService.init();
  }

  public static void main(String[] args) {
    SpringApplication.run(FileuploadApplication.class, args);
  }
  

}

最后启动后端，测试下文件上传。

之前使用的是3D-DNA流程做Hi-C的辅助组装，它的最大优势就是输出结果可以对接下游的JBAT(juicerbox with Assembly Tools)进行手动矫正。然而它点缺陷也很明显，处理速度不够快，且对植物的优化不行，同时目前许久不更新了。

最近我发现了一套新的组合，chromap + yahs 完全替代之前3D-DNA流程。它的依赖工具如下

• chromap: 高效Hi-C数据比对
• samtools: sam转bam
• yahs: 另一个Hi-C scaffolding工具。纠错上准确性高，排序上略强3d-dna，远超SALSA2。
• juicer_tools: 用于输出导入JuiceBox

chrompa, samtools, yahs可以直接用conda进行安装，juicer_tools依赖Java环境，并需要单独下载

1
2
3

conda create -n hic-scaffolding -c bioconda -c conda-forge chromap samtools yahs assembly-stats openjdk 
# 1.19.02版本就行了, 最新的3.0不向下兼容
wget https://s3.amazonaws.com/hicfiles.tc4ga.com/public/juicer/juicer_tools_1.19.02.jar

具体分析步骤如下，我们需要提供前期组装结果，以及Hi-C数据

1
2
3

contigsFasta=/到/你的/contig.fa的路径
r1Reads=/到/你的/Hi-C R1测序的路径
r2Reads=/到/你的/Hi-C R2测序的路径

第一步，数据比对

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 建立索引
samtools faidx $contigsFasta
chromap -i -r $contigsFasta -o contigs.index

# alignment
chromap \
    --preset hic \
    -r $contigsFasta \
    -x contigs.index \
    --remove-pcr-duplicates \
    -1 $r1Reads \
    -2 $r2Reads \
    --SAM \
    -o aligned.sam \
    -t 50

# 排序   
samtools view -bh aligned.sam | samtools sort -@ 50 -n > aligned.bam
rm aligned.sam    

按照read的名字进行排序和按照位置排序或未排序的结果会有一些不同。

第二步，又快又好的scaffolding。默认只需要两个输入，组转的contig.fa和比对的bam，和C语言一样简洁。

1

yahs $contigsFasta aligned.bam

在输出结果中

• inital_break 表示纠错的中间输出
• _scaffolds_final.agp和_scaffolds_final.fa则是最终结果

对于输出结果，我们希望进行可视化，此时可以使用yahs提供的jucier工具

第三步，为juicer_tools准备输入

1
2
3
4
5

juicer pre -a -o out_JBAT \
    yahs.out.bin \
    yahs.out_scaffolds_final.agp \
    $contigsFasta.fai
# -o out_JBAT表示输出文件名的前缀    

一共包括如下几个文件

• out_JBAT.assembly
• out_JBAT.assembly.agp
• out_JBAT.hic
• out_JBAT.liftover.agp
• out_JBAT.txt

out_JBAT.txt就作为下游的输入

1
2
3
4

JUICER=/路径/到/juicer_tools_1.19.02.jar
asm_size=$(awk '{s+=$2} END{print s}' $contigsFasta.fai)
java -Xmx36G -jar $JUICER \
    pre out_JBAT.txt out_JBAT.hic <(echo "assembly ${asm_size}")

输出的out_JBAT.hic就可以导入到JBAT进行组装,导出为out_JBAT.review.assembly

将手动修改的结果传递给juicer，进行scaffolding。

1

juicer post -o out_JBAT out_JBAT.review.assembly out_JBAT.liftover.agp contigs.fa

输出结果为 out_JBAT.FINAL.agp, out_JBAT.FINAL.fa

由于Bioconductor的网站在国外，因此一部分用户可能在使用BiocManager安装R包时因网络问题而失败。

第一个问题是，BiocManager::install()在运行时会下载 http://bioconductor.org/config.yaml 然后检查当前系统是不是符合需求，R包有没有过期。

尽管config.yaml文件非常小，但是依旧有一小部分人会因为访问不了biconductor.org, 导致这步花费非常久的时间。

两个解决方案：

第一个是设置BIOCONDUCTOR_ONLINE_VERSION_DIAGNOSIS来避免自检

1

options(BIOCONDUCTOR_ONLINE_VERSION_DIAGNOSIS=FALSE)

第二个是，下载好这个config.yaml到本地, 设置BIOCONDUCTOR_CONFIG_FILE为该文件在本地的路径

1

options(BIOCONDUCTOR_CONFIG_FILE="/路径/到/你下载的/config.yaml")

第二个问题是，用官方源下载速度非常慢，所以可以使用国内镜像，如

1

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")

亦或者用 chooseBioCmirror() 选择。

如果你拥有最高权限，如果你只管理一台服务器，那么系统自带的包管理工具就帮你解决了所有问题。但是真实世界没有那么美好，所以我花了那么久时间去学习如何从源码开始编译一个软件。

环境为CentOS Linux release 7.4.1708 (Core), Linux内核version 3.10.0-693.el7.x86_64， GCC版本为4.8.5 20150623 (Red Hat 4.8.5-16) (GCC)，

Linux的编译体系

无管理员权限编译的常规三部曲是./configure --prefix=$HOME/usr && make && make install，其中最重要的一步就是configure，它所做的任务如下

• 检查GCC版本以及是否安装了编译所需工具
• 如果需要头文件，则默认去/usr/include查找
• 如果涉及到动态编译库，则默认去/usr/lib和/usr/lib64查找. 注：lib的函数库仅用于开机时用,提供给/bin和/sbin.

那为何需要配置？配置主要解决软件开发和软件实际安装时平台不同所导致的问题。

首先，一个C/C++工程不可能只用到标准库，很多已有的轮子就不需要重复制造。其次，由于很多软件都重复用到相同的依赖库，那么如果把这些依赖库独立成单独的模块，在调用的时候加载，也能节省空间。早期为了适配多个平台，开发人员需要手写条件语句来检查环境依赖，后来GNU专门开发了Autotools辅助构建源码编译所需要的关键文件。

编译环境变量

用./configure -h查看帮助文档的时候，会提示几个和编译相关非常重要的环境变量。

1
2
3
4
5
6
7
8
9
10
11
12

# 编译器
CC          指定C编译器(compiler command)路径
CXX         指定C++编译器
# 编译器选项
CFLAGS      用于C编译器的选项
CXXFLAGS    用于C++编译器的选项
LDFLAGS     链接相关选项,如果你有自定义的函数库(lib dir)，即可以用 -L<lib dir>指定
# 预编译器
CXXCPP      C++ 预处理器(preprocessor)
CPP         C 预处理器(preprocessor)
# 预编译器选项
CPPFLAGS    C/C++预处理器选项, 如果你自定义的头文件，可以用-I<include dir>

Makfile规则中的编译命令通常遵守如下规范：

1,首先从源代码生成目标文件( 预处理,编译,汇编 )，”-c”选项表示不执行链接步骤；

1

$(CC) $(CPPFLAGS) $(CFLAGS) example.c   -c   -o example.o

2,然后将目标文件链接为最终的结果( 链接 )，”-o”选项用于指定输出文件的名字。

1

$(CC) $(LDFLAGS) example.o   -o example

这些只是约定俗成的惯例，所以有些人会“随性而为”，你也拿他没有办法。尽管将源代码编译为二进制文件的四个步骤由不同的程序(cpp,gcc/g++,as,ld)完成，但是事实上 cpp, as, ld 都是由 gcc/g++ 进行间接调用的。换句话说，控制了 gcc/g++ 就等于控制了所有四个步骤。从 Makefile 规则中的编译命令可以看出，编译工具的行为全靠 CC/CXX CPPFLAGS CFLAGS/CXXFLAGS LDFLAGS 这几个变量在控制。所以控制这些变量最简单的做法是首先设置与这些 Makefile 变量同名的环境变量并将它们 export 为 环境变量（全局），然后运行 configure 脚本，大多数 configure 脚本会使用这同名的环境变量代替 Makefile 中的值

• CC/CXX: 指定C/C++编译所在路径，即可以选择不同的版本的编译器进行编译。
• CPPFLAGS: 这是用于预处理阶段的选项。用于添加不在标准路径/usr/include下的头文件。如CPPFLAGS="-I$HOME/usr/include -I$HOME/usr/include/ncurses"
• CFLAGS/CXXFLAGS： 指定头文件（.h文件）的路径，如：CFLAGS=-I/usr/include -I/path/include。同样地，安装一个包时会在安装路径下建立一个include目录，当安装过程中出现问题时，试着把以前安装的包的include目录加入到该变量中来。

CPPFLAG和CFLAGS/CXXFLAGS这两个变量可以认为是等价关系，都用在预处理阶段，也就是都能用于指定头文件所在位置。

• LDFLAGS：gcc 等编译器会用到的一些优化参数，也可以在里面指定库文件的位置。用法：LDFLAGS=-L/usr/lib -L/path/to/your/lib。每安装一个包都几乎一定的会在安装目录里建立一个lib目录。如果明明安装了某个包，而安装另一个包时，它愣是说找不到，可以抒那个包的lib路径加入的LDFALGS中试一下。

有时候LDFLAGS指定-L虽然能让链接器找到库进行链接，但是运行时链接器却找不到这个库，如果要让软件运行时库文件的路径也得到扩展，那么我们需要增加这两个库给”-Wl,R”：

1

LDFLAGS = -L/var/xxx/lib -L/opt/mysql/lib -Wl,R/var/xxx/lib -Wl,R/opt/mysql/lib

如在执行./configure以前设置环境变量 export LDFLAGS="-L/var/xxx/lib -L/opt/mysql/lib -Wl,R/var/xxx/lib -Wl,R/opt/mysql/lib"，注意设置环境变量等号两边不可以有空格，而且要加上引号（shell的用法）。那么执行configure以后，Makefile将会设置这个选项，链接时会有这个参数，编译出来的可执行程序的库文件搜索路径就得到扩展了

除了通过以上几种环境变量为编译器提供头文件和静态和动态库文件的位置信息外，还有一种变量叫做 PKG_CONFIG_PATH , 它从xx.pc文件获取读取相应的环境环境。

注意:Linux下编译共享库时，必须加上-fPIC参数，即export CFLAGS=" -fPIC" CXXFLAGS=" -fPIC"否则在链接时会有错误提示.这是在编译zsh时候发现明明装了ncurse却还是不能用的共享库的坑。

fPIC的目的是什么？共享对象可能会被不同的进程加载到不同的位置上，如果共享对象中的指令使用了绝对地址、外部模块地址，那么在共享对象被加载时就必须根据相关模块的加载位置对这个地址做调整，也就是修改这些地址，让它在对应进程中能正确访问，而被修改到的段就不能实现多进程共享一份物理内存，它们在每个进程中都必须有一份物理内存的拷贝。fPIC指令就是为了让使用到同一个共享对象的多个进程能尽可能多的共享物理内存，它背后把那些涉及到绝对地址、外部模块地址访问的地方都抽离出来，保证代码段的内容可以多进程相同，实现共享。

动态库路径问题: 由前面可以知道许多大型软件为了减少体积不会完全编译所有功能，而是借助于动态连接技术在运行时加载封装好的动态连接库内的功能。这就涉及一个非常重要的问题，软件如何知道动态链接库所处的位置。动态库搜索路径分两种情况，一种是编译生成可执行文件时，另外一种是运行可执行文件时。

编译生成可执行文件时，动态库的搜索路径顺序如下：

• 首先gcc会找-L选项；
• 然后再找gcc的环境变量LIBRARY_PATH，可以在.profile设置这个环境变量，并且可以通过选项-v查看gcc最终编译时LIBRARY_PATH的值；
• 再找内定目录: /lib：/usr/lib：/usr/local/lib，这些都是当初compile gcc时写在程序内的。

注意上面索顺序是不会递归在目录下搜索的。

生成可执行文件后，运行文件时，动态库的搜索路径顺序如下：

• 首先编译目标代码时指定的动态库搜索路径，就是用选项 -Wl,rpath 指定程序在运行时动态库的搜索路径，比如gcc -Wl,-rpath,include -L. -ldltest hello.c，在执行文件时会搜索路径./include；
• 环境变量LD_LIBRARY_PATH指定的动态库搜索路径；
• 配置文件/etc/ld.so.conf中指定的动态库搜索路径，即在配置文件中添加动态库的绝对路径，然后运行指令ldconfig是配置文件生效；
• 默认的动态库搜索路径/lib:/usr/lib。

同样上面索顺序是不会递归在目录下搜索的。通常使用动态库简单做法是：把生成的so文件拷贝到/usr/lib中，这样不管是生成可以执行文件时，还是执行程序时，都能找到需要的so文件。但是普通用户没有/usr/lib的写入权限，所有要指定LD_LIBRARY_PATH.ls

参考资料:

• CFLAGS详解
• Makefile编译选项CC与CXX/CPPFLAGS、CFLAGS与CXXFLAGS/LDFLAGS
• 使用gcc时头文件路径和动态链接库路径

GCC安装(非必要)

首先让我们利用系统原来老旧的GCC编译器编译出最新版本的gcc吧，毕竟安装软件的时候，GCC的版本一定要过最低要求。

第一步： 下载gcc源码

1
2
3
4

mkdir -p ~/src && cd ~/src
wget https://mirrors.tuna.tsinghua.edu.cn/gnu/gcc/gcc-7.2.0/gcc-7.2.0.tar.gz
tar -zxvf gcc-7.2.0.tar.gz && cd gcc-7.2.0
ls

第二步， 检查系统是否已经具备前置软件, 主要是GMP，MPFR, MPC。这些软件可以到ftp://gcc.gnu.org/pub/gcc/infrastructure/找到，然后下载后解压缩，并移动到gcc源码文件夹下。 可以在配置的时候用--with-gmp, --with-mpfr --with-mpc指定具体所在路径。

1
2
3
4
5
6
7
8
9
10
11
12
13

cd src
# GNU Multiple precision Library
wget ftp://gcc.gnu.org/pub/gcc/infrastructure/gmp-6.1.0.tar.bz2 \
&& tar -jxvf gmp-6.1.0.tar.bz2 && mv gmp-6.1.0 gcc-7.2.0/gmp
# isl library
wget ftp://gcc.gnu.org/pub/gcc/infrastructure/isl-0.18.tar.bz2 \
&& tar -jxvf isl-0.18.tar.bz2 && mv isl-0.18 gcc-7.2.0/isl
# MPFR Library
wget ftp://gcc.gnu.org/pub/gcc/infrastructure/mpfr-3.1.4.tar.bz2 \
&& tar -jxvf mpfr-3.1.4.tar.bz2 && mv mpfr-3.1.4 gcc-7.2.0/mpfr
# MPC Library
wget ftp://gcc.gnu.org/pub/gcc/infrastructure/mpc-1.0.3.tar.gz \
&& tar -zxvf mpc-1.0.3.tar.gz && mv mpc-1.0.3 gcc-7.2.0/mpc

不过更加人性化的方法是在GCC源码根目录下运行./contrib/download_prerequisites，可以自动搞定。

第三步：使用./configure进行配置。官方强烈建议, 配置所在文件夹一定要和源码所在文件夹区分开，此外configure还可以配置很多参数，我的代码如下：

1
2
3
4
5

mkdir build && cd build
../configure\
	--prefix=$HOME/usr \ # 指定安装路径
	--disable-multilib \ # 取消32位库编译
	--enable-threads=posix \ # 使用POSIX/Unix98作为线程支持库

基本上这一步不会出现太多的报错，都能够顺利生成Makefile.

第四步： 编译. 这步有一个小技巧就是利用多核处理器进行加速，例如make -j2 就是双核。

这一部分很慢很慢，因为默认设置下是3个阶段的引导(3-stage bootstrap), 以保证能够编译出完整的GCC系统并且还不会出错，你可以在配置的时候用--disable-bootstrap进行关闭。

第五步： 安装。如果你编译都成功了，那么安装也不会存在问题了， make install.

那么我们编译的GCC和系统自带的有什么区别吗？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 从头编译
$ $HOME/usr/bin/gcc -v
Using built-in specs.
COLLECT_GCC=/home/zgxu/usr/bin/gcc
COLLECT_LTO_WRAPPER=/home/zgxu/usr/libexec/gcc/x86_64-pc-linux-gnu/7.2.0/lto-wrapper
Target: x86_64-pc-linux-gnu
Configured with: ../configure --prefix=/home/zgxu/usr --disable-multilib --enable-threads=posix
Thread model: posix
gcc version 7.2.0 (GCC)
# 系统自带
$ gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/libexec/gcc/x86_64-redhat-linux/4.8.5/lto-wrapper
Target: x86_64-redhat-linux
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man --infodir=/usr/share/info --with-bugurl=http://bugzilla.redhat.com/bugzilla --enable-bootstrap --enable-shared --enable-threads=posix --enable-checking=release --with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions --enable-gnu-unique-object --enable-linker-build-id --with-linker-hash-style=gnu --enable-languages=c,c++,objc,obj-c++,java,fortran,ada,go,lto --enable-plugin --enable-initfini-array --disable-libgcj --with-isl=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/isl-install --with-cloog=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-x86_64-redhat-linux/cloog-install --enable-gnu-indirect-function --with-tune=generic --with-arch_32=x86-64 --build=x86_64-redhat-linux
Thread model: posix
gcc version 4.8.5 20150623 (Red Hat 4.8.5-16) (GCC)

不谈安装路径和版本，基本上 差别 就是在配置这一步，而这些参数就需要仔细研究了。

一个 错误 : ‘Link tests are not allowed after GCC_NO_EXECUTABLES.’ 后来发现是第三步没有在独立的文件下构建Makefile.

参考资料：

• installing GCC: https://gcc.gnu.org/install/
• linux下编译gcc6.2.0: https://www.cnblogs.com/oloroso/p/5984985.html

CMake: 平台无关的编译软件

不同平台有着不同的Make工具用于编译,例如 GNU Make ，QT 的 qmake ，微软的 MS nmake，BSD Make（pmake），Makepp，等等。这些 Make 工具遵循着不同的规范和标准，所执行的 Makefile 格式也千差万别。这样就带来了一个严峻的问题：如果软件想跨平台，必须要保证能够在不同平台编译。而如果使用上面的 Make 工具，就得为每一种标准写一次 Makefile ，这将是一件让人抓狂的工作。

CMake就是针对上面问题所设计的工具：它首先允许开发者编写一种平台无关的 CMakeList.txt 文件来定制整个编译流程，然后再根据目标用户的平台进一步生成所需的本地化 Makefile 和工程文件，如 Unix 的 Makefile 或 Windows 的 Visual Studio 工程。从而做到“Write once, run everywhere”。显然，CMake 是一个比上述几种 make 更高级的编译配置工具。一些使用 CMake 作为项目架构系统的知名开源项目有 VTK、ITK、KDE、OpenCV、OSG 等.

1
2
3

wget https://cmake.org/files/v3.10/cmake-3.10.2.tar.gz
tar xf cmake-3.10.2.tar.gz
cd cmake-3.10.2

参考资料：

• http://www.hahack.com/codes/cmake/
• https://www.cnblogs.com/d-blog/p/4617208.html

几个必须要装的函数库

在安装之前需要先声明几个环境变量，可以直接添加在配置文件中。这都是后面遇到共享库的问题得到的经验教训。

1
2
3
4
5
6

export CFLAGS=" -fPIC"
export CXXFLAGS=" -fPIC"
export CPPFLAGS="-I$HOME/usr/include -I$HOME/usr/include/ncurses -I$HOME/usr/include/X11"
export LDFLAGS="-L$HOME/usr/lib -L$HOME/usr/lib64"
export LD_LIBRARY_PATH=$HOME/usr/lib:$HOME/usr/lib64
export PKG_CONFIG_PATH=$HOME/usr/lib/pkgconfig:$HOME/usr/share/pkgconfig

ncurses提供了一系列的函数以便使用者调用它们去生成基于文本的用户界面，许多大名鼎鼎的软件都用到了ncurses，例如vim, screen,tmux,zsh等。并且samtools如果需要tview可视化BAM文件，也需要这个库做支持。

1
2
3

wget ftp://ftp.invisible-island.net/ncurses/ncurses.tar.gz && tar -zxvf ncurses.tar.gz
./configure --enable-shared --prefix=$HOME/usr
make && make install

Libevent是一个用C语言编写的、轻量级的开源高性能事件通知库, 后续安装tmux时候需要这个依赖库。

1
2
3
4
5

# libevent
cd src
wget https://github.com/libevent/libevent/releases/download/release-2.1.8-stable/libevent-2.1.8-stable.tar.gz
tar -zxvf libevent-2.1.8-stable.tar.gz && cd  libevent-2.1.8
./configure prefix=$HOME/usr && make && make install

bzip2, xz, zlib: 文件压缩相关函数库，后续samtools编译时需要。

1
2
3
4
5
6

wget http://www.zlib.net/zlib-1.2.11.tar.gz
tar -zxvf zlib-1.2.11.tar.gz && cd zlib-1.2.11 && ./configure --prefix=$HOME/usr && make && make install
wget http://www.bzip.org/1.0.6/bzip2-1.0.6.tar.gz
tar -zxvf bzip2-1.0.6.tar.gz && cd bzip2-1.0.6 && ./configure --prefix=$HOME/usr && make && make install
wget https://tukaani.org/xz/xz-5.2.3.tar.gz
tar -zxvf xz-5.2.3.tar.gz && cd xz-5.2.3 && ./configure --prefix=$HOME/usr && make && make install

openssl, libssh2, libcurl: 计算机之间文件传输访问相关库。其中OpenSSL是一个安全套接字层密码库，囊括主要的密码算法、常用的密钥和证书封装管理功能及SSL协议，并提供丰富的应用程序供测试或其它目的使用。libssh2是一个C 函数库，用来实现SSH2协议。libcurl主要功能就是用不同的协议连接和沟通不同的服务器.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 安装有先后
# openssl
wget https://www.openssl.org/source/openssl-1.0.2m.tar.gz
tar -zxvf openssl-1.0.2m.tar.gz && cd openssl-1.0.2m
# 这里非常神奇的居然是config，添加shared生成动态库
./config prefix=$HOME/usr shared
make && make install
# 卸载使用 make clean
# libssh2
wget https://www.libssh2.org/download/libssh2-1.8.0.tar.gz
tar -zxvf libssh2-1.8.0.tar.gz && cd libssh2-1.8.0
./configure --with-libssl-prefix=$HOME/usr/ssl --prefix=$HOME/usr
# libcurl
wget https://curl.haxx.se/download/curl-7.56.1.tar.gz
tar -zxvf curl-7.56.1.tar.gz && cd curl-7.56.1
./configure --prefix=$HOME/usr --enable-http --enable-ftp --enable-file --enable-proxy --enable-telnet --enable-libcurl-option --enable-ipv6 --with-lib --with-ssl

readline: GNU提供用于这些命令补全、搜索历史命令、行编辑快捷键等等这些人性化的交互方式的函数库，缺少这个标准库，编译的R就缺少自动补全的功能。

1
2
3

wget http://ftp.gnu.org/gnu/readline/readline-7.0.tar.gz
tar -zxvf readline-7.0.tar.gz && cd readline-7.0
./configure --prefix=$HOME/usr && make && make install

PCRE: 提供和Perl5相同语法和语义正则表达式的函数库，后续安装R用到。

1
2
3

wget ftp://ftp.csx.cam.ac.uk/pub/software/programming/pcre/pcre-8.41.tar.gz
tar -zxvf pcre-8.41.tar.gz && cd pcre-8.41
./configure --enable-utf --enable-pcregrep-libz --enable-pcregrep-libbz2 --prefix=$HOME/usr

x11：X11也叫做X Window系统，X Window系统 (X11或X)是一种位图显示的视窗系统,是在 Unix 和 类Unix 操作系统，以及OpenVMS上建立图形用户界面的标准工具包和协议，并可用于几乎所有已有的现代操作系统。主要是R编译的时候要用，具体用途有待开发。

x11安装比较复杂，有很多的依赖库，因此需要先安装xtrans, xextproto, xcb(lib,xcb-proto, libpthread-subs), kbproto,xproto,inputproto。但是编译很容易，仅提供下载链接

1
2
3
4
5
6
7
8
9
10
11

https://www.x.org/releases/X11R7.7/src/lib/xtrans-1.2.7.tar.gz
https://www.x.org/releases/X11R7.7/src/proto/xextproto-7.2.1.tar.gz
https://www.x.org/releases/X11R7.7/src/proto/kbproto-1.0.6.tar.gz
https://www.x.org/releases/X11R7.7/src/proto/xproto-7.0.23.tar.gz
https://www.x.org/releases/X11R7.7/src/proto/inputproto-2.2.tar.gz
https://www.x.org/releases/X11R7.7/src/xcb/libpthread-stubs-0.3.tar.gz
https://www.x.org/releases/X11R7.7/src/xcb/xcb-proto-1.7.1.tar.gz
https://www.x.org/releases/X11R7.7/src/xcb/libxcb-1.8.1.tar.gz
https://www.x.org/releases/X11R7.7/src/lib/libSM-1.2.1.tar.gz
https://www.x.org/releases/X11R7.7/src/lib/libICE-1.0.8.tar.gz
https://www.x.org/releases/X11R7.7/src/lib/libXt-1.1.3.tar.gz

相当于人工检查依赖环境，仅仅是繁琐而已，并不困难

1
2
3
4

# 安装X11
wget -4 https://www.x.org/releases/X11R7.7/src/lib/libX11-1.5.0.tar.gz
tar -zxvf libX11-1.5.0.tar.gz && cd libX11-1.5.0
./configure --prefix=$HOME/usr && make && make install

编译案例

安装zsh

zsh或许可以认为是最好的shell，用过zsh的人都不会想着bash了。不过zsh自定义配置，可供选择的插件以及主题实在是太多，因此一定要搭配oh-my-zsh。zsh依赖ncurses.

1
2
3
4
5

wget -O zsh.tar.gz https://sourceforge.net/projects/zsh/files/latest/download
tar -zxvf zsh.tar.gz && cd zsh
export CPPFLAGS="-I$HOME/usr/include/" LDFLAGS="-L$HOME/usr/lib"
../configure --prefix=$HOME/usr --enable-shared
make && make install

由于没有root权限，无法使用chsh，只能通过在~/.bashrc添加exec $HOME/usr/bin/zsh -l保证登陆的时候自动切换成zsh。其次, zsh搭配oh-my-zsh才完整, 只不过这里只能手动安装了。

1
2
3
4
5
6
7
8

# 从github上克隆oh-my-zsh
git clone git://github.com/robbyrussell/oh-my-zsh.git ~/.oh-my-zsh
# 用oh-my-zsh的zsh配置文件替代
cp ~/.oh-my-zsh/templates/zshrc.zsh-template ~/.zhsrc
# 安装一些字体, 不然一些主题会显示异常
cd src
git clone https://github.com/powerline/fonts.git --depth=1
cd fonts && ./install.sh

重启一下终端，后面根据需要调整配置文件里的参数。

编译tmux

tmux和screen类似，也是文本终端神器, 依赖于libevent和ncurses.

1
2
3
4
5
6
7
8

export CPPFLAGS="-I$HOME/usr/include -I$HOME/usr/include/ncurses"
export LDFLAGS="-L$HOME/usr/lib -L$HOME/usr/lib64"
mkdir -p src && cd src
git clone https://github.com/tmux/tmux.git
cd tmux
sh autogen.sh
./configure --prefix=$HOME/usr
make && make install

编译R语言

由于我自己编译完全版的GCC套餐，很多之前的gfortran不存在的问题也就不存在了（管理员安装了Java）。此外，R还需要gnu readline, pcre > 8.2, x11。当然这些函数包都在之前安装好了。

一些依赖库

1
2

# 安装tidyverse发现xm12需要libiconv的libiconv.so.2
https://ftp.gnu.org/pub/gnu/libiconv/libiconv-1.15.tar.gz

正式安装

1
2
3
4

wget https://cran.r-project.org/src/base/R-3/R-3.4.2.tar.gz
tar -zxvf R-3.4.2.tar.gz  && cd R-3.4.2/
./configure --prefix=$HOME/R
make && make install

到此，我可以说Linux平台下即便我没有root权限，也没有多少软件包是我所不能手工编译。

最新上线chatGPT，非常的热闹，openAI的服务器几度因为请求太多而宕机，紧急增加了服务器。很多视频都说程序员岌岌可危，因为chatGPT真的能帮我们写代码。

我试了试，发现chatGPT确实好用，比如说我让他帮我写一段差异分析的代码(differential analysis in R code)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# Load the necessary packages
library(DESeq2)

# Import the RNA-seq data into R
data <- readRDS("rnaseq_data.rds")

# Create a DESeqDataSet object from the data
dds <- DESeqDataSet(data, design = ~ condition)

# Perform the differential expression analysis
dds <- DESeq(dds)

# Extract the results of the analysis
results <- results(dds)

# Filter the results to only include significantly differentially expressed genes
results <- results[which(results$padj < 0.05), ]

这是他的实现方式，确实有模有样。但其实，我如果去查DESeq2的相关的教程，也能找到类似的代码呀。

另外，有些时候，它还会给出错误的实现方式，比如说，我让他绘制GWAS的Q-Q图，以下是它给的代码。

1
2
3
4
5
6
7
8
9

# Load the necessary packages
library(graphics)
library(GWASdata)

# Load the example GWAS data
data(gwas1)

# Create a Q-Q plot for the GWAS data
qqplot(gwas1$pvalue)

其中有一个R包，GWASdata需要额外安装。于是，我问他怎么安装呢？他给出的建议是用 install.packages("GWASdata")。然而，事实上，GWASdata是放在Bioconductor上的，所以安装方式是BiocManager::install("GWASdata")。当然即便我安装上了，但是这个代码也无法运行，因为gwas1这个数据集不存在与GWASdata中，我目前也不知到和这个gwas1哪里来的。更有趣的是，经过我的检索，实际上，我们应该用qqPlot这个来自于qqman包的函数，才能根据p值来画图。

实际上，用AI辅佐写代码的事情，之前用GitHub的Copilot时，已经感受到了它的实用性了。它能够帮忙补全很多代码，然后我们在此基础上做些修改，差不多就能用了。但是，真的让它完成整个项目，我觉得还有一段路要走吧。

事情源于，我在写脚本的时候，在脚本里面插入了一句 conda activate 环境名, 然后出现如下的报错提示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

    $ conda init <SHELL_NAME>

Currently supported shells are:
  - bash
  - fish
  - tcsh
  - xonsh
  - zsh
  - powershell

See 'conda init --help' for more information and options.

IMPORTANT: You may need to close and restart your shell after running 'conda init'.

conda认为我没有初始化环境，我脚本是在bash环境下运行的，我也用conda init bash初始化过。因此，问题肯定不是出在这里。

利用关键词”conda activate in bash script”检索，我找到了两种解决方法。

方法1: 在脚本中多加一句

1

source $HOME/miniconda/etc/profile.d/conda.sh

注意，我的conda是安装在家目录下的miniconda目录中，对于非家目录的安装方式，要修改 $HOME/miniconda。

方法2: 我们可以通过 conda run 来运行给定环境下的命令，假如，我们安装了一个环境rna-seq, 里面有一个程序叫做STAR, 我们可以随便写一个tmp.sh脚本，内容为

1

conda run -n rna-seq STAR --help

那么，此时运行 bash tmp.sh 就不会报错。也就是说，你并不是一定要用conda activate 启动环境，才能调用命令，你其实可以调用某个环境的给定指令。

方法2相对于方法1有个非常大的优势，那就是，如果你有多个不同python版本的环境，你不用担心写脚本的时候写了启动，但是忘了写退出。你只需要在原来的代码前加上一句， conda run -n 环境名。

• calling conda source activate from bash script
• activating conda environment in within a shell script