抢救你破碎的测序数据

生信技能树

发布于 2022-06-27 20:45:22

4450

发布于 2022-06-27 20:45:22

文章被收录于专栏：生信技能树生信技能树

最近一个优秀学员写了笔记给我，所以就多聊了几句，关心一下对方学生信的目的以及现有的ngs项目，结果就挖掘到了一个难题：

我下载了对方传递给我的如下所示的文件：

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校验了，本来是想把它们全部先解压再说，马上就报错：


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命令：

$ 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文件，如下所示：

$ 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，如下所示：

# 第一个文件：
$ 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，其实也足够定量。所以我写了如下所示的取巧的代码：

 ls *fq |while read id; do (head -n 30000000 $id|gzip > ${id}.gz );done

这样，不管每个样品的文库大小如何，也不管它是何种程度的破碎，我猜测它起码前面的 7.5M的测序reads是ok的，所以就对它们进行下面的定量流程。得到如下所示的文件：

$ 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文件：

$ 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文件，后面就是比对到参考基因组，并且对基因进行定量即可，如下所示：

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的转录组测序的推荐数据量，但是做差异分析理论上也足够啦。后续差异分析非常简单了，公众号推文在：

我简单跑了一下，确实没有问题，一个简单的火山图，如下所示：

火山图

因为测序数据量确实不够，所以我们的流程里面过滤低表达量基因后就只剩下1.2万个基因左右啦，如果是标准的20M的转录组测序的推荐数据量，火山图里面通常是有2~3万个基因，甚至加大测序量还可以探索编码和非编码。甚至融合基因事件，可变剪切。不过现在我们就抢救到了少量数据，仅仅是能大致保证差异分析是问题不大。

但是，这个抢救你破碎的测序数据过程其实需要两个前提：