新冠肺炎CT辅助诊断文献实战-01

Good

新冠肺炎CT辅助诊断文献实战-01

   https://www.nature.com/articles/s41551-020-00633-5
2020年11月，华中科技大学发表在Nature Biomedical Engineering
如果不想看文献的话可以看我写的文献导读（文献导读001）基于图像识别技术的cov19辅助诊断 - 知乎 (zhihu.com)
建议还是看一下，文章不难，我也只用了30分钟就大概的看了一遍
  由于上一个系列——深度学习与单细胞，我自己在写教程运行的过程中发现了一系列的问题导致目前暂时无法更新。然后最近老师推荐看来一篇关于医学图像识别辅助诊断新冠肺炎的文章，看了一下文章发现特别简单，于是就打算利用这篇文章继续我的深度学习系列教程。
   这个系列教程与深度学习和生物信息学的应用暂时无关，就是非常基础的pytorch实战和一系列的问题。
网络上其实可以搜索到很多相关的教程，只是我这个是针对医学数据而已哈
  数据下载

   iCTCF - CT images and clinical features for COVID-19 (biocuckoo.cn)
  一共有两套数据需要下载：

• 已经标记好的label数据，用于构建第一个模型
  
• cohort1 和cohort2
  

cohort1，cohort2的下载略微有点麻烦得一个个点，这里的话可以使用迅雷自带的批量下载功能
如果是ubuntu系统的话，可以使用如下代码批量下载:

1. # 这里本身就是国内源，不要做奇奇怪怪的操作
2. for i in {1..1521}  
3. do  
4.      wget -c http://ictcf.biocuckoo.cn/patient/CT/Patient%20{$i}.zip
5. done

复制代码

背景

根据文章的神经网络框架，该框架分为三个部分其中，第一部分和第二部分是用图像数据训练，第三部分是用临床特征数据进行训练。
由于CT图像一次性会扫描出很多张图片，但是并不是每一张图片都是有用的。因此，医生在对CT图像做出判断决策时候通常需要以下几个步骤：

• 从大量的CT图像中挑选出包涵信息的图片
  
• 从这些图片中挑选出可能与疾病相关的图片
  
• 从可能与疾病相关的图片中做出决策
  

上面的三个步骤又可以总结成为两个步骤，即：

• 挑选图片
  
• 做出决策
  

那么深度学习技术的辅助诊断也是从这两个部分进行的。那么，在完全不考虑经费和计算资源的情况下那么辅助诊断的实施就非常的简单。直接整合所有的CT图像数据和clinical outcomes 建模即可，但是，显然这种做法并不现实，因为CT图像通常一个样本产生的数据非常多，并且由于噪声越大，所需要的模型也越大，训练时长也越久。因此，文章的做法是将两个模型分开做，即构建一个挑选图片的模型，将挑选后的图片再次构建一个clinical outcomes 相关的模型，这个构思也是目前辅助诊断中最为常用的模型。此外，将步骤分开的话，由于单个任务的复杂度低，因此构建模型也相对简单。
正文

那么，我们先按照文献给的思路一步步的构建模型，先构建第一个模型
   文章使用的是kears，我用的是pytorch哈，由于文章中有些部分我感觉并不是很合理，因此我打算稍微改动一下
  首先初始化配置

1. #### no mean but to fold the code
2. import os
3. import math
4. import itertools
5. import numpy as np
6. import pandas as pd2
7. import matplotlib.pyplot as plt
8. from collections import Counter
9. from matplotlib import cm
10. import cv2
12. import torch
13. import torchvision
14. import torch.utils.data as Data
15. import torch.nn as nn
16. import torch.nn.functional as F
17. import pytorch_lightning as pl
18. from sklearn.preprocessing import StandardScaler,MinMaxScaler
19. from sklearn.model_selection import train_test_split
21. import warnings
22. from myexptorch import *
23. %matplotlib inline
24. %config InlineBackend.figure_format = 'svg'
25. warnings.filterwarnings("ignore")
27. # if GPU avaliabel, data to
28. device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
29. def deviceChange(x):
30.     if (torch.cuda.is_available() & x == 'cuda'):
31.         device = torch.device('cuda')
32.     else:
33.         device = torch.device('cpu')
34.     return device
36. # =============================================================================
37. # the following code is common pytorch formula
38. # =============================================================================
40. # change expr to tensor
41. def expToTorch(x):
42.     return torch.from_numpy(np.array(x)).float()
45. # change label to tensor
46. def labelTotorch(y):
47.     return torch.LongTensor(y)
49. # data to data iter
50. def makeDataiter(x,y,batch_size,shuffle=True):
51.     return Data.DataLoader(Data.TensorDataset(x, y), batch_size, shuffle=shuffle)
53. # scale
54. def DataFrame_normalize(my_data,std_method):
55.     if std_method == 'minmax':
56.         method = MinMaxScaler()
57.     if std_method == 'std':
58.         method = StandardScaler()
59.     my_data = pd.DataFrame(method.fit(my_data).transform(
60.         my_data),columns=my_data.columns,index=my_data.index)
61.     return my_data
63. def toOneHot(ylabels,n_class):
64.     onehot = torch.zeros(ylabels.shape[0],n_class)
65.     index = torch.LongTensor(ylabels).view(-1,1)
66.     onehot.scatter_(dim=1, index=index, value=1)
67.     return onehot
69. def toLabel(ylabels):
70.     return torch.topk(ylabels, 1)[1].squeeze(1)
72. # show tenor to picture
73. def show_tensor_pictures(x, y):
74.     plt.figure(figsize=(12,12))
75.     fig, axes = plt.subplots(1, len(x))
76.     for ax, image, label in zip(axes, x, y):
77.         ax.imshow(image.view((28, 28)).numpy())
78.         ax.set_title(label)
79.         ax.axes.get_xaxis().set_visible(False)
80.         ax.axes.get_yaxis().set_visible(False)
81.     plt.show()
83. def show_batch(imgs):
84.     grid = torchvision.utils.make_grid(imgs,nrow=5)
85.     plt.imshow(grid.numpy().transpose((1, 2, 0)))
86.     plt.title('Batch from dataloader')

复制代码

导入数据

pytorch的数据导入方法非常的方便，还会自动添加label，导入的数据仅需要满足如下要求：

• 创建一个用于储存数据的文件夹 (totallabeldata)
  
• 在该文件夹下将同label的数据存放在相同的文件夹下，并命名
  

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-x33KBSOd-161********44)(https://tva1.sinaimg.cn/large/008eGmZEly1gnylnfyvh2j310r06ggm9.jpg)]
然后使用如下代码进行数据导入

1. # load images to torch
2. # plz put same class images to dif-folder
3. train = torchvision.datasets.ImageFolder('../totallabeldata/',
4.                                             transform=torchvision.transforms.Compose([
5.                                                 torchvision.transforms.Grayscale(num_output_channels=1),
6.                                                 torchvision.transforms.Scale(256),
7.                                                 torchvision.transforms.CenterCrop(200),
8.                                                 torchvision.transforms.ToTensor()])
9.                                             )
10. train_loader = torch.utils.data.DataLoader(train, batch_size=128,shuffle=True,num_workers =24)
11. len(train_loader.dataset.targets)
12. train_loader.dataset.class_to_idx
14. """
15. 在train_loader.dataset.class_to_idx的对象中记录了每个类别对应的class number
16. """
17. # [out] :
18. # 19685
19. # {'NiCT': 0, 'nCT': 1, 'pCT': 2}

复制代码

然后这里的话就面临着一个问题，就是没办法进行data的分割，总所周知，我们在进行建模的时候，通常要将数据分割成3份：training, validation, test，通常情况下比例为8:1:1。 因此要对数据进行如下操作

1. # 目前我没有找到更好的方法了，如果有知道的朋友请留个言呗
2. # Create empty list
3. XRawData = []
4. YRawData = []
6. # from iteration extract batch
7. for i,(x,y) in enumerate(train_loader):
8.     XRawData.append(x)
9.     YRawData.append(y)
11. # concat all batch
12. XRawData = torch.cat(XRawData)
13. YRawData = torch.cat(YRawData)
15. randoms = 42
16. batch_size =150
17. # Step.3 data split to train - test
18. x_train, x_test, y_train, y_test = train_test_split(
19.     XRawData, YRawData, stratify = YRawData,
20.     test_size=0.1, random_state=randoms)
22. # Step.4 train data split validation data
23. x_train, x_val, y_train, y_val = train_test_split(
24.     x_train, y_train, stratify = y_train,
25.     test_size=0.1, random_state=randoms)
27. # Step.6 transform to iterate object
28. train_iter = makeDataiter(x_train, y_train, batch_size=batch_size,shuffle=True)
29. val_iter = makeDataiter(x_val, y_val, batch_size=batch_size,shuffle=True)
31. # visualization
32. for i, (batch_x, batch_y) in enumerate(train_iter):
33.     if(i<1):
34.         show_batch(batch_x)
35.         plt.axis('off')
36.         plt.show()
38. # class count
39. Counter(YRawData.numpy())
40. # [out]:
41. # Counter({2: 4001, 1: 9979, 0: 5705})

复制代码

构建模型（按文章）

1. class VGG_Simple(pl.LightningModule):
3.     train_epoch_loss = []
4.     train_epoch_acc = []
5.     train_epoch_aucroc = []
7.     val_epoch_loss = []
8.     val_epoch_acc = []
9.     val_epoch_aucroc = []
11.     test_predict = []
12.     test_sample_label = []
13.     test_decoder = []
15.     test_digitcaps = []
16.     test_matrix = []
17.     test_conv = []
18.     test_primary = []
19.     test_decoder = []
21.     def __init__(self,myloss = None):
22.         super(VGG_Simple, self).__init__()
23.         self.myloss = myloss
24.         self.cov = nn.Sequential(
25.             nn.Conv2d(in_channels=1,out_channels=64,kernel_size=3,stride=1),
26.             nn.ReLU(inplace=True),
27.             nn.Conv2d(in_channels=64,out_channels=64,kernel_size=3,stride=1),
28.             nn.ReLU(inplace=True),
29.             nn.MaxPool2d((2, 2)),
31.             nn.Conv2d(in_channels=64,out_channels=32,kernel_size=3,stride=1),
32.             nn.ReLU(inplace=True),
33.             nn.Conv2d(in_channels=32,out_channels=32,kernel_size=3,stride=1),
34.             nn.ReLU(inplace=True),
35.             nn.MaxPool2d((2, 2)),   
37.             nn.Conv2d(in_channels=32,out_channels=16,kernel_size=3,stride=1),
38.             nn.ReLU(inplace=True),
39.             nn.Conv2d(in_channels=16,out_channels=16,kernel_size=3,stride=1),
40.             nn.ReLU(inplace=True),
41.             nn.MaxPool2d((2, 2))      
42.         )
44.         self.dnn = nn.Sequential(
45.             nn.Linear(7056, 64),
46.             nn.ReLU(inplace=True),
47.             nn.Dropout(0.5),
48.             nn.Linear(64, 32),
49.             nn.ReLU(inplace=True),
50.             nn.Dropout(0.5),
51.             nn.Linear(32, 3)
52.         )
54.     def forward(self, x):
56.         # x.size(0)即为batch_size
57.         in_size = x.size(0)
58.         out = self.cov(x)
59.         out = out.view(in_size, -1)
60.         out = self.dnn(out)
61.         return out
63.     # 3. 定义优化器
64.     def configure_optimizers(self):
65.         optimizer = torch.optim.Adam(self.parameters(), lr=0.001,weight_decay=0.05)
66.         return optimizer
68.     # 4. 训练loop
69.     def training_step(self, train_batch, batch_ix):
70.         vector_sample, sample_label = train_batch
71.         probs = self.forward(vector_sample)
72.         loss  = self.myloss(probs, sample_label)
73.         acc = pl.metrics.functional.accuracy(probs, sample_label)
74.         mylogdict = {'loss': loss, 'log': {'train_loss': loss, 'train_acc': acc}}
75.         return mylogdict
77.     def validation_step(self, validation_batch, batch_ix):
78.         vector_sample, sample_label = validation_batch
79.         probs= self.forward(vector_sample)
80.         loss  = self.myloss(probs, sample_label)
81.         val_acc = pl.metrics.functional.accuracy(probs, sample_label)
82.         self.log_dict({'val_loss': loss, 'val_acc': val_acc})
83.         mylogdict = {'log': {'val_loss': loss, 'val_acc': val_acc}}
84.         return mylogdict
86.     def test_step(self, test_batch, batch_ix):
87.         vector_sample, sample_label = test_batch
88.         probs= self.forward(vector_sample)
89.         self.test_predict.append(probs.cpu())
90.         self.test_sample_label.append(sample_label.cpu())
91.         return {'test': 0}
93.     def training_epoch_end(self, output):
94.         train_loss = sum([out['log']['train_loss'].item() for out in output]) / len(output)
95.         self.train_epoch_loss.append(train_loss)
97.         train_acc = sum([out['log']['train_acc'].item() for out in output]) / len(output)
98.         self.train_epoch_acc.append(train_acc)
99.         print(train_acc)
100.         return train_loss
102.     def validation_epoch_end(self, output):
103.         val_loss = sum([out['log']['val_loss'].item() for out in output]) / len(output)
104.         self.val_epoch_loss.append(val_loss)
106.         val_acc = sum([out['log']['val_acc'].item() for out in output]) / len(output)
107.         self.val_epoch_acc.append(val_acc)
108.         print('mean_val_loss: ', val_loss, '\t', 'mean_val_acc: ', val_acc)
109.         return val_loss

复制代码

然后开始训练模型

1. # training
2. traning = True
3. # set seed
4. pl.utilities.seed.seed_everything(seed=42)
6. # Defined function
7. myloss = nn.CrossEntropyLoss()
8. model = VGG_Simple(myloss=myloss)
10. if traning:
11.     # output file
12.     OUTPUT_DIR = './lightning_logs'
13.     tb_logger = pl.loggers.TensorBoardLogger(save_dir='./',
14.                                              name=f'check_point')
16.     # set check point to choose best model
17.     checkpoint_callback = pl.callbacks.ModelCheckpoint(
18.         dirpath=tb_logger.log_dir,
19.         filename='{epoch}-{val_acc:.5f}-{val_loss:.5f}',
20.         save_top_k = 15, #保留最优的15
21.         monitor='val_acc', # check acc
22.         mode='auto'
23.     )
25.     # train loop
26.     trainer = pl.Trainer(gpus=-1,
27.                          callbacks=[checkpoint_callback],
28.                          max_epochs = 100,
29.                          gradient_clip_val=0,
30.                          auto_lr_find=False)# 开始训练
31.     trainer.fit(model, train_iter, val_iter)
33.     from pathlib import Path
34.     out = Path(tb_logger.log_dir)
35.     print(out)
36.     [ckpt.stem for ckpt in out.iterdir()]

复制代码

但是，如果你这样训练那么会出现如下的问题：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-p21vXGnK-161********46)(https://tva1.sinaimg.cn/large/008eGmZEly1gnylnhvxatj30u00t642p.jpg)]
不管如何训练，loss不会改变，accuracy也不会改变，那么这是什么原因呢？常见的有如下两种原因：

• 模型复杂度不够
  
• 模型一开始就过拟合
  

当然，在这次训练中，这两种情况可以完全排除，理由很简单，我用的模型和文章一模一样
那么，我们就要从loss函数或者数据本生出发，这里我们要注意一下小细节，在之前我们进行class counter的时候的输出

1. # class count
2. Counter(YRawData.numpy())
3. # [out]:
4. # Counter({2: 4001, 1: 9979, 0: 5705})

复制代码

发现，样本量的差距着实有点大，但是说大把也不是很大，可能在运行三五百次的epoch应该是完全没有什么太大的问题
但是为了节约时间我们要对loss进行加权

1. class_weight = torch.FloatTensor([2,5,2]).cuda()
2. myloss = nn.CrossEntropyLoss(weight=class_weight)

复制代码

最终修改代码如下

1. # training traning = True# set seed pl.utilities.seed.seed_everything(seed=42)# Defined function class_weight = torch.FloatTensor([2,5,2]).cuda()
2. myloss = nn.CrossEntropyLoss(weight=class_weight)
3. model = VGG_Simple(myloss=myloss)if traning:    # output file    OUTPUT_DIR = './lightning_logs'    tb_logger = pl.loggers.TensorBoardLogger(save_dir='./',                                             name=f'check_point')    # set check point to choose best model     checkpoint_callback = pl.callbacks.ModelCheckpoint(        dirpath=tb_logger.log_dir,        filename='{epoch}-{val_acc:.5f}-{val_loss:.5f}',        save_top_k = 15, #保留最优的15        monitor='val_acc', # check acc         mode='auto'    )    # train loop     trainer = pl.Trainer(gpus=-1,                          callbacks=[checkpoint_callback],                         max_epochs = 100,                         gradient_clip_val=0,                         auto_lr_find=False)# 开始训练    trainer.fit(model, train_iter, val_iter)    from pathlib import Path    out = Path(tb_logger.log_dir)    print(out)    [ckpt.stem for ckpt in out.iterdir()]

复制代码

输出结果如下

我们可以看到，在第8次epoch就已经达到92%的准确率了，但是，还有一个问题，这个准确率还是太低了。
考虑到，CT数据的图像是灰度图，且图像的大致形状也差不多，感觉各个医院拍出来的图也差不了多少（自己认为的），所以我认为，模型中的dropout层是没有用的
   这里是我觉得文章不妥的第一个点

1. model.train_epoch_acc[40:44]
2. # [out]
3. # [ 0.8633503256557143,
4. # 0.8609685636012354,
5. # 0.862724444576513,
6. # 0.8625120386899074]

复制代码

我们发现train的过程准确率也才86%
在删除dropout后，准确率得到显著的提升

1. model.train_epoch_acc[40:44]
2. # [out]
3. #[0.9764486033225728,
4. # 0.9797252973663473,
5. # 0.979190********92,
6. # 0.9742679161446117]

复制代码

验证模型

1. # load model parameter
2. version = 'version_14' # choose best model path
3. path = os.path.join(os.getcwd(),'check_point',version)
4. ckpt = str(os.listdir(path)[2]) # choose best model
5. print(ckpt)
7. # get best model parameter
8. ckpt = os.path.join(path,ckpt)
9. model = VGG_Simple(myloss=myloss)
10. trainer = pl.Trainer(resume_from_checkpoint=ckpt, gpus=-1)
12. # create test_iter
13. test_iter = makeDataiter(x_test, y_test, batch_size=batch_size,shuffle=False)
14. trainer.test(model, test_iter)
16. # get test accuracy
17. predict = torch.cat(model.test_predict)
18. y = torch.cat(model.test_sample_label)
19. predict = toLabel(predict).cpu().numpy()
20. y = y.cpu().numpy()
21. sum( y == predict) / len(y)
23. # [out]:
24. # 0.9862874555611986

复制代码

测试集准确率98%，然后这里要注意一下，文章用的是AUC去评估，对于像这种的3分类任务，使用AUC评估很不直观，应该使用混淆矩阵或者accuarcy进行评估
   这里是我觉得文章不妥的第二个点

1. Best Regards,  
2. Yuan.SH;
3. School of Basic Medical Sciences,  
4. Fujian Medical University,  
5. Fuzhou, Fujian, China.  
6. please contact with me via the following ways:  
7. (a) e-mail :yuansh3354@163.com  

复制代码

来源：https://blog.csdn.net/qq_40966210/article/details/114021421
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！