view release on metacpan or  search on metacpan

doc/masurca.md  view on Meta::CPAN

├── translateReduceFile.perl
└── ufasta

0 directories, 100 files
```

同时还生成一个配置文件样例, `sr_config_example.txt`.

# 样例数据

MaSuRCA 发表在 Bioinformatics 时自带的测试数据.

> IMPORTANT! Do not pre‐process Illumina data before providing it to MaSuRCA. Do not do any
> trimming, cleaning or error correction. This WILL deteriorate the assembly

Super-reads在 `work1/superReadSequences.fasta`, `work2/` 和 `work2.1/` 是 short jump 的处理, 不用管.
`superReadSequences_shr.frg` 里面的 super-reads 是作过截断处理的, 数量不对.

> Assembly result. The final assembly files are under CA/10-gapclose and named 'genome.ctg.fasta'
> for the contig sequences and 'genome.scf.fasta' for the scaffold sequences.

MaSuRCA-3.1.3 supports gzipped fastq files while MaSuRCA-2.1.0 doesn't.

## Rhodobacter sphaeroides (球形红细菌)

高 GC 原核生物 (68%), 基因组 4.5 Mbp.

```bash
mkdir -p ~/data/test
cd ~/data/test

wget -m ftp://ftp.genome.umd.edu/pub/MaSuRCA/test_data/rhodobacter .

mv ftp.genome.umd.edu/pub/MaSuRCA/test_data/rhodobacter .
rm -fr ftp.genome.umd.edu
find . -name ".listing" | xargs rm
```

### Illumina PE, Short Jump and Sanger (1x or 4x)

```bash
cd ~/data/test

cat <<EOF > sr_config.txt
PARAMETERS
CA_PARAMETERS = ovlMerSize=30 cgwErrorRate=0.25 merylMemory=8192 ovlMemory=4GB 
LIMIT_JUMP_COVERAGE = 60
KMER_COUNT_THRESHOLD = 1
EXTEND_JUMP_READS = 0
NUM_THREADS = 16
JF_SIZE = 50000000
END

EOF

# Illumina PE, Short Jump and Sanger4
mkdir -p rhodobacter_PE_SJ_Sanger4
cp sr_config.txt rhodobacter_PE_SJ_Sanger4/
cat <<EOF >> rhodobacter_PE_SJ_Sanger4/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
JUMP= sj 3600 200  /home/wangq/data/test/rhodobacter/SJ/short_1.fastq  /home/wangq/data/test/rhodobacter/SJ/short_2.fastq
OTHER=/home/wangq/data/test/rhodobacter/Sanger/rhodobacter_sphaeroides_2_4_1.4x.frg
END

EOF

# Illumina PE, Short Jump and Sanger
mkdir -p rhodobacter_PE_SJ_Sanger
cp sr_config.txt rhodobacter_PE_SJ_Sanger/
cat <<EOF >> rhodobacter_PE_SJ_Sanger/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
JUMP= sj 3600 200  /home/wangq/data/test/rhodobacter/SJ/short_1.fastq  /home/wangq/data/test/rhodobacter/SJ/short_2.fastq
OTHER=/home/wangq/data/test/rhodobacter/Sanger/rhodobacter_sphaeroides_2_4_1.1x.frg
END

EOF

# Illumina PE and Short Jump
mkdir -p rhodobacter_PE_SJ
cp sr_config.txt rhodobacter_PE_SJ/
cat <<EOF >> rhodobacter_PE_SJ/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
JUMP= sj 3600 200  /home/wangq/data/test/rhodobacter/SJ/short_1.fastq  /home/wangq/data/test/rhodobacter/SJ/short_2.fastq
END

EOF

# Illumina PE, and Sanger4
mkdir -p rhodobacter_PE_Sanger4
cp sr_config.txt rhodobacter_PE_Sanger4/
cat <<EOF >> rhodobacter_PE_Sanger4/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
OTHER=/home/wangq/data/test/rhodobacter/Sanger/rhodobacter_sphaeroides_2_4_1.4x.frg
END

EOF

# Illumina PE, and Sanger
mkdir -p rhodobacter_PE_Sanger
cp sr_config.txt rhodobacter_PE_Sanger/
cat <<EOF >> rhodobacter_PE_Sanger/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
OTHER=/home/wangq/data/test/rhodobacter/Sanger/rhodobacter_sphaeroides_2_4_1.1x.frg
END

EOF

# Illumina PE
mkdir -p rhodobacter_PE_Sanger
cp sr_config.txt rhodobacter_PE_Sanger/
cat <<EOF >> rhodobacter_PE_Sanger/sr_config.txt
DATA
PE=  pe 180 20 /home/wangq/data/test/rhodobacter/PE/frag_1.fastq /home/wangq/data/test/rhodobacter/PE/frag_2.fastq
END

EOF

# Run
cd ~/data/test

for d in rhodobacter_PE_SJ_Sanger4 rhodobacter_PE_SJ_Sanger rhodobacter_PE_SJ rhodobacter_PE_Sanger4 rhodobacter_PE_Sanger rhodobacter_PE rhodobacter_superreads;
do
    echo "==> ${d}"
    if [ -e ${d}/work1/superReadSequences.fasta ];
    then
        continue     
    fi

    pushd ~/data/test/rhodobacter_PE_SJ_Sanger4 > /dev/null
    $HOME/share/MaSuRCA/bin/masurca sr_config.txt
    bash assemble.sh
    popd > /dev/null
done
```

### Rhodobacter sphaeroides with `anchr superreads`

```bash
# gzip original fastq
mkdir -p ~/data/test/rhodobacter/PEgz
gzip -c ~/data/test/rhodobacter/PE/frag_1.fastq > ~/data/test/rhodobacter/PEgz/frag_1.fq.gz
gzip -c ~/data/test/rhodobacter/PE/frag_2.fastq > ~/data/test/rhodobacter/PEgz/frag_2.fq.gz

mkdir -p ~/data/test/rhodobacter_superreads
cd ~/data/test/rhodobacter_superreads

perl ~/Scripts/sra/superreads.pl \
    ~/data/test/rhodobacter/PEgz/frag_1.fq.gz \
    ~/data/test/rhodobacter/PEgz/frag_2.fq.gz \
    -s 180 -d 20

```

### 结果比较

```bash
cd ~/data/test/

printf "| %s | %s | %s | %s | %s | %s | %s | %s |\n" \
    "Name" "N50SR" "#SR" "N50Contig" "#Contig" "N50Scaffold" "#Scaffold" "EstG" \
    > stat.md
printf "|:--|--:|--:|--:|--:|--:|--:|--:|\n" >> stat.md

for d in rhodobacter_PE_SJ_Sanger4 rhodobacter_PE_SJ_Sanger rhodobacter_PE_SJ rhodobacter_PE_Sanger4 rhodobacter_PE_Sanger rhodobacter_PE rhodobacter_superreads;
do
    printf "| %s | %s | %s | %s | %s | %s | %s | %s |\n" \
        ${d} \
        $( faops n50 -H -N 50 -C ${d}/work1/superReadSequences.fasta ) \
        $( faops n50 -H -N 50 -C ${d}/CA/10-gapclose/genome.ctg.fasta ) \
        $( faops n50 -H -N 50 -C ${d}/CA/10-gapclose/genome.scf.fasta ) \
        $( cat ${d}/environment.sh \
            | perl -n -e '/ESTIMATED_GENOME_SIZE=\"(\d+)\"/ and print $1' )
done >> stat.md

cat stat.md
```

| name          | N50SR |  #SR | N50Contig | #Contig | N50Scaffold | #Scaffold |    EstG |
|:--------------|------:|-----:|----------:|--------:|------------:|----------:|--------:|
| PE_SJ_Sanger4 |  4586 | 4187 |    205225 |      69 |     3196849 |        35 | 4602968 |
| PE_SJ_Sanger  |  4586 | 4187 |     63274 |     141 |     3070846 |        28 | 4602968 |
| PE_SJ         |  4586 | 4187 |     43125 |     219 |     3058404 |        59 | 4602968 |
| PE_Sanger4    |  4705 | 4042 |    125228 |      67 |      534852 |        30 | 4595684 |
| PE_Sanger     |  4705 | 4042 |     19435 |     412 |       21957 |       359 | 4595684 |
| PE            |  4705 | 4043 |     20826 |     407 |       34421 |       278 | 4595684 |
| superreads    |  4705 | 4043 |           |         |             |           | 4595684 |

有足够多的 long reads 支持下, 不需要 short jump.

# SuperReads 3.1.3

2017 年 2 月, UMD ftp 上多了一个新程序
[SuperReads_RNA](ftp://ftp.genome.umd.edu/pub/MaSuRCA/beta/SuperReads_RNA-1.0.1.tar.gz), 是 MaSuRCA
3.2.1 的简化版. 很可能是 `StringTie` 用了 super-reads 来处理 RNA-seq, 在很多人的要求下做的.

根据这个版本, 我将 MaSuRCA 3.1.3 简化, 去掉所有的依赖, 去掉配合 `Celera Assembler` 的部分, 只留下了
`SuperReads`, 可以用 `Linuxbrew` 安装.

```bash
brew install homebrew/science/jellyfish
brew install wang-q/tap/quorum@1.1.1
brew install wang-q/tap/superreads
```

# Super-reads and anchors

## E. coli: link anchors

```bash