最近一个优秀学员写了笔记给我，所以就多聊了几句，关心一下对方学生信的目的以及现有的ngs项目，结果就挖掘到了一个难题：
我下载了对方传递给我的如下所示的文件：
代码语言：javascript复制ls -lh *gz|cut -d" " -f 5-
2.2G 5月 27 19:50 KO_4_1.fq.gz
2.2G 5月 27 19:50 KO_4_2.fq.gz
2.8G 5月 27 19:50 KO_5_1.fq.gz
2.4G 5月 27 19:51 KO_5_2.fq.gz
2.2G 5月 27 19:51 KO_6_1.fq.gz
2.1G 5月 27 19:51 KO_6_2.fq.gz
2.1G 5月 27 19:51 WT_1_1.fq.gz
1.5G 5月 27 19:51 WT_1_2.fq.gz
2.5G 5月 27 19:51 WT_2_1.fq.gz
2.0G 5月 27 19:51 WT_2_2.fq.gz
2.7G 5月 27 19:51 WT_3_1.fq.gz
2.1G 5月 27 19:51 WT_3_2.fq.gz
肉眼看起来没有啥问题，因为对方本来就是测序数据文件破损了，所以也无所谓md5校验了，本来是想把它们全部先解压再说，马上就报错：
代码语言：javascript复制
ls *gz |xargs gunzip
gzip: KO_4_1.fq.gz: invalid compressed data--crc error
gzip: KO_4_1.fq.gz: invalid compressed data--length error
既然 gunzip 命令不支持这样的破碎的测序数据文件，我只能是退而求其次，使用zcat命令：
代码语言：javascript复制$ ls *gz|while read id;do(zcat $id > ${id%%.*}.fq );done
gzip: KO_4_1.fq.gz: invalid compressed data--crc error
gzip: KO_4_1.fq.gz: invalid compressed data--length error
它这个时候虽然也会报错，但是好歹是可以输出内容，得到的全部的fq文件，如下所示：
代码语言：javascript复制$ ls -lh *fq|cut -d" " -f 5-
4.9G 5月 27 19:54 KO_4_1.fq
4.3G 5月 27 19:54 KO_4_2.fq
4.5G 5月 27 19:55 KO_5_1.fq
5.3G 5月 27 19:56 KO_5_2.fq
4.9G 5月 27 19:56 KO_6_1.fq
4.1G 5月 27 19:57 KO_6_2.fq
5.1G 5月 27 19:57 WT_1_1.fq
3.8G 5月 27 19:58 WT_1_2.fq
3.1G 5月 27 19:58 WT_2_1.fq
4.2G 5月 27 19:59 WT_2_2.fq
3.0G 5月 27 19:59 WT_3_1.fq
4.3G 5月 27 20:00 WT_3_2.fq
$ wc -l *fq
57077430 KO_4_1.fq
50696111 KO_4_2.fq
52266180 KO_5_1.fq
61837264 KO_5_2.fq
57638234 KO_6_1.fq
47515910 KO_6_2.fq
59866115 WT_1_1.fq
44587122 WT_1_2.fq
35635542 WT_2_1.fq
49281737 WT_2_2.fq
35259629 WT_3_1.fq
50382180 WT_3_2.fq
可以看到，虽然每个样品都是双端测序，所以都是两个fq文件，但是都不完整！我们通过取巧的方式，可以解压出来一些测序的reads，如下所示：
代码语言：javascript复制# 第一个文件：
$ head KO_4_1.fq
@ST-E00192:686:HL5FYCCXY:7:1101:11048:1538 1:N:0:NGTACTAG
CCTTTCCTGCATGTCACAGCAGAAAAACACAGTGCTCGGAGGT
AAA-AA<FFJFJFFJFJJ7F<FJ<JJJFJJJA<-F-F-77F<7<
@ST-E00192:686:HL5FYCCXY:7:1101:13261:1538 1:N:0:NGTACTAG
TCCCAGATTTCTGTAGCAGATGATGATGATGAAAGTCTTCTAGGACATT
AAFFAA7FAFJJA-<-A<AF-7F<<FF-FF-FA-7-A--<<FF-<<A---
# 第二个文件：
$ head KO_4_2.fq
@ST-E00192:686:HL5FYCCXY:7:1101:11048:1538 2:N:0:NGTACTAG
NGCTACGACGTGCTCATCCAGAAGTTCCTGAGCCTGTACGGC
#-A<FJJJJAAAF<JFJJJFJFAAA<JJ<J<-FF<FF7A<AJ
@ST-E00192:686:HL5FYCCXY:7:1101:13261:1538 2:N:0:NGTACTAG
NCTCCATATACAGGCAAAACTACGTGTGCCATTTATCTTAC
#AAFFJFFJJJJAFFFFJJJFJFAF7--AJAJA<JFJ-<FF
现在的问题就是，同一个样品的双端测序数据的两个fq文件里面需要保证reads的一致性。
所以，首先需要遍历同一个样品的双端测序数据的两个fq文件，拿到那些匹配的ID，然后按照顺序输出成为两个fq文件，这样它们的reads数量就相等，而且顺序还是一致的。但是我已经七八年没有写脚本了，不管是perl还是Python，所以这里就不献丑了，非常简单其实，可以作为一个生信工程师考核题。
不过，凭借我的聪明才智，我这里猜测到了另外一个取巧的办法，其实我们的转录组测序数据量都很大，20M的reads绰绰有余，而我们前面的 wc -l *fq 发现大家都是大于30M的行，而30M的行起码还有 7.5M的测序reads，因为一个测序reads会占用4行。而 7.5M的测序reads，其实也足够定量。所以我写了如下所示的取巧的代码：
代码语言：javascript复制 ls *fq |while read id; do (head -n 30000000 $id|gzip > ${id}.gz );done
这样，不管每个样品的文库大小如何，也不管它是何种程度的破碎，我猜测它起码前面的 7.5M的测序reads是ok的，所以就对它们进行下面的定量流程。得到如下所示的文件：
代码语言：javascript复制$ ls -lh *gz|cut -d" " -f 5-
570M 5月 27 20:13 KO_4_1.fq.gz
591M 5月 27 20:17 KO_4_2.fq.gz
689M 5月 27 20:21 KO_5_1.fq.gz
585M 5月 27 20:24 KO_5_2.fq.gz
602M 5月 27 20:28 KO_6_1.fq.gz
578M 5月 27 20:31 KO_6_2.fq.gz
566M 5月 27 20:35 WT_1_1.fq.gz
633M 5月 27 20:38 WT_1_2.fq.gz
749M 5月 27 20:43 WT_2_1.fq.gz
583M 5月 27 20:46 WT_2_2.fq.gz
759M 5月 27 20:51 WT_3_1.fq.gz
569M 5月 27 20:54 WT_3_2.fq.gz
这个时候跑trim_galore确实也没有问题，得到如下所示的过滤后的fq文件：
代码语言：javascript复制$ ls -lh 2.clean_fq/*val*|cut -d" " -f 5-
552M 5月 27 21:12 2.clean_fq/KO_4_1_val_1.fq.gz
569M 5月 27 21:12 2.clean_fq/KO_4_2_val_2.fq.gz
660M 5月 27 21:13 2.clean_fq/KO_5_1_val_1.fq.gz
561M 5月 27 21:13 2.clean_fq/KO_5_2_val_2.fq.gz
583M 5月 27 21:12 2.clean_fq/KO_6_1_val_1.fq.gz
558M 5月 27 21:12 2.clean_fq/KO_6_2_val_2.fq.gz
547M 5月 27 21:12 2.clean_fq/WT_1_1_val_1.fq.gz
608M 5月 27 21:12 2.clean_fq/WT_1_2_val_2.fq.gz
721M 5月 27 21:14 2.clean_fq/WT_2_1_val_1.fq.gz
561M 5月 27 21:14 2.clean_fq/WT_2_2_val_2.fq.gz
725M 5月 27 21:14 2.clean_fq/WT_3_1_val_1.fq.gz
547M 5月 27 21:14 2.clean_fq/WT_3_2_val_2.fq.gz
可以看到，同样的测序数据，同一个样品过滤前后，其实变化并不大，主要是因为测序已经是比较稳定的技术啦。
有了过滤好的fq文件，后面就 是比对到参考基因组，并且对基因进行定量即可，如下所示：
代码语言：javascript复制Sample Name % Assigned M Assigned
KO_4.sort 69.2% 5.9M
KO_5.sort 63.0% 5.5M
KO_6.sort 69.7% 6.0M
WT_1.sort 67.7% 5.8M
WT_2.sort 73.3% 6.1M
WT_3.sort 68.4% 5.9M
因为每个样品都是 7.5M的测序reads，所以最后的定量也是在6M附近，它虽然达不到20M的转录组测序的推荐数据量，但是做差异分析理论上也足够啦。后续差异分析非常简单了，公众号推文在：
- 解读GEO数据存放规律及下载，一文就够
- 解读SRA数据库规律一文就够
- 从GEO数据库下载得到表达矩阵 一文就够
- GSEA分析一文就够（单机版 R语言版）
- 根据分组信息做差异分析- 这个一文不够的
- 差异分析得到的结果注释一文就够
我简单跑了一下， 确实没有问题，一个简单的火山图，如下所示：
火山图
因为测序数据量确实不够，所以我们的流程里面过滤低表达量基因后就只剩下1.2万个基因左右啦，如果是标准的20M的转录组测序的推荐数据量，火山图里面通常是有2~3万个基因，甚至加大测序量还可以探索编码和非编码。甚至融合基因事件， 可变剪切。不过现在我们就抢救到了少量数据，仅仅是能大致保证差异分析是问题不大。
但是，这个抢救你破碎的测序数据过程其实需要两个前提：
- 首先你破碎的不能太严重
- 其次破碎的发生是随机的，但是不破坏reads顺序
