Viral-Track： 揭示重症COVID-19病人的感染特征
前言
今年的新冠可谓是来得触不及防，对各个行业多带来了很大的影响。庆幸的是在全国人民一起配合的努力下，国内的疫情被控制的很好，相比于海外来说，简直是天壤之别了。现在想想能够作为华夏人民的一员还真的是很荣幸！
既然疫情来了，那么就勇敢面对吧！在面对的同时，大家也肯定想弄清楚疫情爆发的原因。对于一探究竟的任务，那是科研工作的事。作为吃瓜群众，我们只需坐等别人公布答案就可以了。今天想跟大家分享的是一篇方法学的文章《Host-Viral Infection Maps Reveal Signatures of Severe COVID-19 Patients》，这篇文章是5月份发表在cell杂志上的，该文章虽然不能告诉你新冠爆发的原因，但却可以用来检测病毒的感染情况来揭示重症病人的感染特征。这其中的方法还是值得借鉴和学习的。值得提醒的是这个方法通用的检测其他病毒，还是挺厉害的。
我并不想跟大家解读文献，感兴趣的可以点击文章链接自己详细阅读，这里主要展示一下分析过程。
准备工作
准备好scRNAseq的数据，就是把read1中的cell barcode和umi序列追加到read2中序列名的后面，然后根据推测的barcode白名单来抽取出reads序列，这样所有信息都在read2中，后续分析就是用这个处理后的read2。完成这一步可以使用umi_tools软件来处理，如果不会使用这款工具的小伙伴可以看我之前分享的帖子《umi_tools for dealing with UMIs and cell barcodes》。不过得提醒一下，这个软件不支持多线程，数据量大的话会很耗时。来直观地看一下处理前后数据的变化，由于后面只用到处理后的read2，我这里就只展示read2的前后变化了。
- 处理前的序列
@SRX8108992.1 1 length=151
NATAAATATCACCAGTAATAAGATGTATGGTATCATAAATGCTTTGATATTATCCATTGAGAAGGATCCACCATATTTTTTGTTGAATTGCCAAAAATGCATAACTTCATTCCAATCATAAAACNCTGGAGAAACCCCAAGGGAAGGACGT
+
#<A<<<JFJJFJJFFJFJJJFAJJFJJJJJFJJJ7J-AJJJJJFJ<7F-FJ-AJJJ<FJ-FJ<JJJJFJAJJ-J-<7FJJJJJ-FA7AJJJJJJJ<AA7--A<-7FAFAAJJJJ7AF<JF77-<#--<<-)--7<<F)--77)--7A<77)
@SRX8108992.2 2 length=151
NTGTTGGTGAAGGAAGAAGTGGGGTGGAAGAAGTGGGGTGGGACGACAGTGAAATCTAGAGTAAAACCAAGCTGGCCCAAGGTGTCCTGCAGGCTGTAATGCAGTTTAATCAGAGTGCCATTTTNTTTTTGTTCCAATGATTTTAATTATT
+
#<A<FJJJJJAJJ<-FFFJJJJJJAFJAAFFFJJJJJJJJJJ<JJJJFJJJFFAFJJFJ7FJ7F<FJFFFFJJJJFJ<A-<F--AFFJJFJJJJ7<JJ<-FF---AJF<FJ7JFFJJJJJJFJF#FFFAAAFA<----<7<A----77<-7
- 处理后的序列
@SRX8108992.2840_TCTGAGAAGAGTACCG_TGGTCCTTCG 2840 length=151
GCAGTTTTTTGAAAATGGGCTCAACCAGAAAAGCCCAAGTTCATGCAGCTGTGGCAGAGTTACAGTTCTGTGGTTTCATGTTAGTTGCCTTATAGTTACAGTGTAATTAGTGCCACTTAATGTATGTTACCAAACATTAATATATATACCC
+
AAAF<FJJ<JAJ<F-FAFAF7FJFJJFJAJJJFJFJJJAJFJJFJJAF<<F-FJJ-F-FAAA<-7FFJ7FAFJF-7A<-FFJ-7AF-AFJJ7F<FAJ-7-AFJ-7<FFAFJFAF<JJJF-FFJFJ<<F---A-77<F-7A-<-<7-<F-AF
@SRX8108992.2841_ATTTCTGAGTATCGAA_CCACACATTG 2841 length=151
ATGAATTCCTTGCCTGATATTGAAGAAGTCAAGGATAAAAGCATAAAGAAAGGAAATAATGCCTTCAGGGTCTACCGAATGCTGCCCCTATCAGAACGGCCTTCTAAGAAAGGAAAGAAACCAAAGACAGAAAAAGACGACAAAGTTAAGA
+
AAA<FFJJJ77FFJJJJJJJJJJJJJJJJJJFAJJJJJJJJJFJ<A<A<JJJJFJJJJJFJFFJ<JFAJJ-FFJAAJJF<JJJFJJJFJFAJJ-<FJFF--7FJ--JJJJJAFFJFJ<FFFA---7-<FAFJF<-7-----<-<F--7A<-
大家注意到变化了么，序列长度没有任何变化，只是序列的名字多了barcode和umi信息，基因过滤掉了那些不是barcode白名单中的序列。到此可用于后续分析的read2文件就准备好了。接下来就是走流程了，这个前提是在你的环境中已经安装好需要的R包("Biostrings","ShortRead","doParallel","GenomicAlignments","Gviz", "GenomicFeatures","Rsubread")和需要用的软件(STAR、samtools、StringTie)。
1、Detection of viruses in scRNA-seq data
第一步是准备病毒和宿主的参考基因组的STAR比对索引文件，病毒基因组可从网站ViruSite上下载，宿主基因组可从 ensembl上下载。构建比对索引命令如下：
STAR --runThreadN N --runMode genomeGenerate --genomeDir /path/to/index --genomeFastaFiles /path/to/Virusite_file.fa path/to/Host_genome_chromosome*.fa
第二步准备病毒的注释文件(制表符分隔)，内容如下：
Name_sequence Genome_length Virus_name Complete_segment_name
NC_007646 135135 Ovine herpesvirus 2 strain BJ1035 Ovine herpesvirus 2 strain BJ1035, complete genome.
NC_003401 133719 Macacine herpesvirus 5 Macacine herpesvirus 5, genome.
NC_008725 111953 Maruca vitrata MNPV Maruca vitrata MNPV, complete genome.
NC_001777 3684 Tobacco necrosis virus A Tobacco necrosis virus A, complete genome.
NC_004203 1119 Banna virus strain JKT-6423 Banna virus strain JKT-6423, segment 8
NC_003521 241087 Panine herpesvirus 2 strain Heberling Panine herpesvirus 2 strain Heberling, complete genome.
NC_016440 8638 Garlic common latent virus Garlic common latent virus, complete genome.
NC_016562 132662 Campylobacter phage CPX Campylobacter phage CPX, complete genome.
第三步准备好两个配置文件，Parameters.txt 和 Files_to_process.txt 分别用来存放脚本需要的配置和待检测的fastq文件路径，内容分别如下：
# Parameters.txt 配置文件内容如下：
N_thread=8
Output_directory="/home/chenyl/project/scRNAseq_COVID-19_result"
Index_genome="/home/chenyl/database/virusite/star_virusite_hg38"
Viral_annotation_file="/home/chenyl/database/virusite/virusite_anno.txt"
Name_run="covid19" #任务名，可自行取名字
Load_STAR_module=FALSE
Load_samtools_module=FALSE
Load_stringtie_module=FALSE
Minimal_read_mapped=50
Single_cell_metadata="" #此参数可用为空
# Files_to_process.txt 配置文件内容如下：
/home/chenyl/fastq/SRX8108992_extracted_2.fastq.gz
/home/chenyl/fastq/SRX8108993_extracted_2.fastq.gz
/home/chenyl/fastq/SRX8108994_extracted_2.fastq.gz
上面的 Parameters.txt配置文件中的参数，大家看到前面的变量名估计就知道是设置什么的了，这里我就不解释了，只需要把变量后面的内容替换为自己的内容就可以了。
最后就是来运行流程脚本了，命令如下：
nohup Rscript $scriptdir/Viral_Track_scanning.R Parameters.txt Files_to_process.txt >log 2>&1 &
这个检测脚本需要运行很长时间(STAR比对不是并发比对的，而是一个样本一个样本的循环比对；还有对bam文件的一些统计)，大家还是放在后台让它自己慢慢运行吧。
等待程序运行完毕，且没有报错的情况下，会在结果目录中给每个样本生成一个存放结果文件的目录(提示一下‘QC_report.pdf’是检测的pdf报告)，每个样本的目录结构类似如下：
SRX8108994_extracted_2
├── Count_chromosomes.txt
├── Merged_viral_mapping.bam
├── QC_filtered.txt
├── QC_report.pdf
├── QC_unfiltered.txt
├── SRX8108994_extracted_2_Aligned.sortedByCoord.out.bam
├── SRX8108994_extracted_2_Aligned.sortedByCoord.out.bam.bai
├── SRX8108994_extracted_2_Log.final.out
├── SRX8108994_extracted_2_Log.out
├── SRX8108994_extracted_2_Log.progress.out
├── SRX8108994_extracted_2_SJ.out.tab
├── SRX8108994_extracted_2__STARgenome
│ ├── sjdbInfo.txt
│ └── sjdbList.out.tab
├── SRX8108994_extracted_2__STARpass1
│ ├── Log.final.out
│ └── SJ.out.tab
└── Viral_BAM_files
├── GL000009.2.bam
├── GL000224.1.bam
├── KI270707.1.bam
├── KI270717.1.bam
├── KI270726.1.bam
├── KI270728.1.bam
├── KI270734.1.bam
├── KI270747.1.bam
├── refseq|NC_001479|7835nt|Encephalomyocarditis.bam
├── refseq|NC_001782|1801nt|Saccharomyces.bam
├── refseq|NC_004102|9646nt|Hepatitis.bam
├── refseq|NC_008310|6485nt|Hibiscus.bam
├── refseq|NC_008580|8380nt|Rabbit.bam
├── refseq|NC_009127|295146nt|Cyprinid.bam
├── refseq|NC_009758|8926nt|Marine.bam
├── refseq|NC_009823|9711nt|Hepatitis.bam
├── refseq|NC_013221|2896nt|Phytophthora.bam
├── refseq|NC_022518|9472nt|Human.bam
├── refseq|NC_028095|1705nt|Torulaspora.bam
├── refseq|NC_029302|6038nt|Piscine.bam
├── refseq|NC_031032|24320nt|Bacillus.bam
├── refseq|NC_038425|9531nt|Non-primate.bam
├── refseq|NC_038828|1866nt|Heterobasidion.bam
├── refseq|NC_040589|11413nt|Diatom.bam
├── refseq|NC_040615|166748nt|Eptesicus.bam
├── refseq|NC_041925|54836nt|Proteus.bam
└── refseq|NC_043054|137452nt|Bubaline.bam
QC_report.pdf报告展示如下：
基于上面的检测是主要的分析内容，下面还有两个可选的分析，这两个分析内容我没有实际跑过分析脚本，大家如果感兴趣可以自己动手分析一下。
2、Transcriptome assembly
可以组装检测的病毒转录本，生成病毒的GTF文件，命令如下：
Rscript Viral_Track_transcript_assembly.R Parameter_file.txt Files_to_process.txt
3、Single-cell demultiplexing
提取病毒在每一个细胞中的基因表达量生成一个表格，可用于一步的分析，命令如下：
Rscript Viral_Track_cell_demultiplexing.R Parameter_file.txt Files_to_process.txt
最后
emm，今天就分享到这里，各位看官们支持一下顺便点个赞吧！！！
