1. SOPRA 简介

SOPRA(Statistical Optimization of Paired Read Assembly)，主要用于利用成对的reads信息来进行scaffold连接，其准确性非产高。其参考文献：Dayarian, Adel, Todd P. Michael, and Anirvan M. Sengupta. “SOPRA: Scaffolding algorithm for paired reads via statistical optimization.” BMC bioinformatics 11.1 (2010): 345.

2. SOPRA的下载和安装

$ wget http://www.physics.rutgers.edu/~anirvans/SOPRA/SOPRA_v1.4.6.zip
$ mkdir /opt/biosoft/SOPRA_v1.4.6/
$ unzip SOPRA_v1.4.6.zip -d /opt/biosoft/SOPRA_v1.4.6/
安装目录/opt/biosoft/SOPRA_v1.4.6/下有SOPRA的manual文档
$ chmod 755 /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
$ perl -p -i -e 's\s*$/\n/' /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
SOPRA的主要运行程序就是这4支perl程序。

3. SOPRA的使用流程

SOPRA 可以支持 Illumina 和 SOLID 的二代测序数据，但此处仅讲解使用 Illumina 的测序数据。SOPRA 也支持同时利用多个paired测序文库的数据，同时支持paired-end和mate-paired数据。

3.1 准备输入数据

不管是paired-end或mate-paired数据，首先需要将数据的方向转换为FR方向。即需要将mate-paired的2端reads都进行反向互补。然后，需要将fatstq文件2端的reads文件合并为一个文件，并转换为fasta文件作为SOPRA的输入数据。
此外，SOPRA还需要contig序列作为输入的基因组序列。
以1个paired-end数据和1个mate-paired数据为例：

先将 jumping 文库的数据进行反向互补
$ fastq_rc.pl jump.1.fastq > 11; mv 11 jump.1.fastq
$ fastq_rc.pl jump.2.fastq > 22; mv 22 jump.2.fastq

然后，合并两端的reads序列
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  frag.1.fastq frag.2.fastq frag.fastq
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  jump.1.fastq jump.1.fastq jump.fastq

再将fastq转换为fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' frag.fastq > frag.fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' jump.fastq > jump.fasta

3.2 SOPRA 对contigs序列和测序数据进行格式处理

使用s_prep_contigAseq_v1.4.6.pl程序对基因组的contig序列以及illunia测序数据进行序列名的修改，以利于后续程序运行。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_prep_contigAseq_v1.4.6.pl \
  -contig contig.fasta -mate frag.fasta jump.fasta -a SOPRA_OUT/

程序参数：
-contig string
    输入基因组 contigs 序列。
-mate string
    输入 illumina 数据。如果有多个数据，则多个fasta文件之间使用空格分开。
-a string
    设置输出文件夹目录。

在 SOPRA_OUT 目录中生成了 contigs_sopra.fasta、frag_sopra.fasta 和 jump_sopra.fasta 文件。查看这3个文件内容，发现是对fasta文件中的序列名进行了修改。

3.3 将Illumina数据比对到基因组

可以利用bowite,bwa等序列比对软件将Illumina数据按单端序列比对到基因组序列上。若一条序列比对到多个位置，推荐报告多个结果。因为后续的s_parse_sam_v1.4.6.pl程序会去除掉比对到多个位置的reads。若仅报告最优结果，我不知道程序会如何处理。要是不去除这样的比对结果，可能会对scaffolding的准确性有影响。因此，还是要注意报告多个比对结果，特别是bowtie2默认下仅报告最优结果。
此外，去除reads尾部碱基质量低的碱基后，illumina reads 的匹配率若是提高很多，则推荐去除reads尾部10~20个碱基，再用于比对。由于mate-pair数据进行了反向互补，则其原本的尾部成了首部。

$ cd SOPRA_OUT/
$ bowtie2-build contigs_sopra.fasta contig
$ bowtie2 -p 20 -x contig -k 10 -f -3 15 -U frag_sopra.fasta -S frag_sopra.sam
$ bowtie2 -p 20 -x contig -k 10 -f -5 15 -U jump_sopra.fasta -S jump_sopra.sam

程序参数：
-k 10
    若read比对到多个位点，则最多报告 10 个结果。
-3 15
    去除read尾部15bp碱基后再用于比对。
-5 15
    去除read首部15bp碱基后再用于比对。

3.4 对SAM文件进行分析

使用s_parse_sam_v1.4.6.pl对sam文件进行分析并简化其信息。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_parse_sam_v1.4.6.pl \
  -sam frag_sopra.sam jump_sopra.sam -a ./

程序参数：
-sam string
    输入 sam 文件。多个illumina 文库的sam文件用空格隔开。
-a string
    设置输出文件的路径。

得到sam文件分析后的结果frag_sopra.sam_parsed和jump_sopra.sam_parsed。

3.5 进行序列方向和距离分析

读取上一步简化后的sam文件信息，进行分析。得到contigs序列的覆盖度、library文库的插入片段长度，序列长度等。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_read_parsed_sam_v1.4.6.pl \
  -parsed frag_sopra.sam_parsed -d 500 -parsed jump_sopra.sam_parsed -d 3000 -a ./

程序参数：
-parsed string
    输入 .sam_parsed 文件。若有多个文件，则用空格分割。
-d int
    illumina数据的插入片段长度。
-c int default: 5
    如果read及其反向互补序列在library中出现的次数>=该值，则不计算其成对信息。
    在运行 s_scaf_v1.4.6.pl 程序后，输入日志中得到如下信息：
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100，则需要考虑降低 -c 参数的值。
    此外，如果数据量比较少，平均覆盖度较低，比如低于 10，则考虑增加 -c 参数的值。

-e int default：0
    如果设置其值为 1， 则使用输入的插入片段长度值。默认下软件会计算插入片段长度值，并使用软件计算出来的值。
-a string
    设置输出文件的路径。

输出文件为 orientdistinfo_c5 。 此外，在 div 文件夹中有contigs序列的覆盖度、library文库的插入片段长度，序列长度等统计信息。

3.6 进行scaffold连接

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_scaf_v1.4.6.pl \
  -o orientdistinfo_c5 -a ./

程序常用参数：
-o string
    输入文件 orientdistinfo_c。
-w int default: 4
    连接2个contigs所需要的最小的pair links。
    在运行 s_scaf_v1.4.6.pl 程序后，输入日志中得到如下信息：
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100，则需要考虑提高 -w 参数的值。一般 -w 取该值的 4~5%。如果该值较低，比如低于10的时候，则需要降低 -w 参数的值。
-L int default: 150
    用于连接scaffold的最小长度的contigs序列。
-h float default: 2.2
    设置一个系数，用于排除对高覆盖度contigs的scaffold连接。大于 mean_coverage + h * std_coverage 该覆盖度的 contigs ，不能用于scaffold连接
-a string
    设置输出文件的路径。

程序输出文件为 scaffolds_h2.2_L150_w4.fasta 。该文件为最终的结果文件。