#!/usr/bin/env bash

## 2017-07-03
# - testing clean-up for the 2-cell dataset
#$> sshi #(interactive node @ cori)

thisDP='/global/projectb/sandbox/rqc/ayspunde/work_sandbox/jira/seqqc-jira/seqqc-13976/GNTNS'
bam_in_FP_2138='/global/dna/dm_archive/sdm/pacbio/00/21/38/pbio-2138.19153.bc1004.bc1004_BAK8A--bc1004_BAK8A.bam'
bam_in_FP_2140='/global/dna/dm_archive/sdm/pacbio/00/21/40/pbio-2140.19198.bc1004.bc1004_BAK8A--bc1004_BAK8A.bam'

thisOutDP="$thisDP"

for this_bam_FP in $bam_in_FP_2138 $bam_in_FP_2138; do
  bam_in_FN=$(basename $this_bam_FP)
  bam_out_w_FP="$thisDP"/"$bam_in_FN"
  bam_out_b_FP=$(echo $bam_out_w_FP | sed 's/\.bam/\.b\.bam/')
  date; time icecreamfinder.sh \
      in=$this_bam_FP out=$bam_out_w_FP outb=$bam_out_b_FP \
      cq=f ow=t keepshortreads=f sambamba=f
  
  samtools sort \
    -o $bam_out_w_FP -O bam -n -T temp_sort -@ 16 \
     $bam_out_w_FP
  $SMRT_BIN_701/pbindex $bam_out_w_FP
done

exit 0

new_ds_name='GNTNS_01'
cd $thisOutDP
xml_2138=$(echo $bam_in_FP_2138 | sed 's/\.bam/\.subreadset\.xml/')
xml_2140=$(echo $bam_in_FP_2140 | sed 's/\.bam/\.subreadset\.xml/')

# #$ Failed: no valid bam files ..
# $SMRT_BIN_701/dataset merge \
#     $new_ds_name.'subreadset.xml' \
#     $xml_2138 $xml_2140

$SMRT_BIN_701/pbindex 'pbio-2138.19153.bc1004.bc1004_BAK8A--bc1004_BAK8A.bam'
$SMRT_BIN_701/pbindex 'pbio-2140.19153.bc1004.bc1004_BAK8A--bc1004_BAK8A.bam'

$SMRT_BIN_701/dataset create \
  --type SubreadSet \
  --name $new_ds_name \
  --generateIndices \
    $new_ds_name.subreadset.xml \
    'clean_input.fofn' 
  
# validate 
$SMRT_BIN_701/dataset validate $new_ds_name.subreadset.xml
## perfecto: PB1205_vanilla_004.subreadset.xml is valid DataSet XML with valid ResourceId references

# import
SMRT_HOST_701='dint01.nersc.gov'  # 'https://dint01.nersc.gov:8243'
SMRT_PORT='9091'
XML2IMP=$thisOutDP/$new_ds_name.subreadset.xml
# CMD
$SMRT_BIN_701/pbservice import-dataset \
  --host $SMRT_HOST_701 --port $SMRT_PORT -u admin -p admin \
"$XML2IMP"

new_ds_name='GNTNS_01M'
$SMRT_BIN_701/dataset create \
  --type SubreadSet \
  --name $new_ds_name \
  --generateIndices \
  --metadata "/global/dna/dm_archive/sdm/pacbio/00/21/38/pbio-2138.19153.bc1001.bc1001_BAK8A--bc1001_BAK8A.subreadset.xml" \
    $new_ds_name.subreadset.xml \
    'clean_input.fofn' 

*DATASET SUMMARY*
          id: 4672
        uuid: 6b9a279f-8467-d554-9e30-3379dcdc1408
        name: PB1205_vanilla_004
  numRecords: 0
 totalLength: 0
   createdAt: 2019-06-25T19:39:36.000Z
   updatedAt: 2019-06-25T19:39:36.000Z
  importedAt: 2019-06-26T02:45:50.962Z
   createdBy: admin
        tags:
    isActive: true
 parent UUID:
 numchildren: 0
       jobId: 3883
   projectId: 1
        path: /global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/PB1205_out_0625/PB1205_vanilla_004.subreadset.xml


*JOB SUMMARY*
            id: 3883
          uuid: f63428d5-32eb-4707-b980-f2d4b6227f04
          name: import-dataset
         state: SUCCESSFUL
    project id: 1
     jobTypeId: import-dataset
     is active: true
     createdAt: 2019-06-26T02:45:49.374Z
     updatedAt: 2019-06-26T02:45:51.001Z
  jobUpdatedAt: 2019-06-26T02:45:51.001Z
      run time: 2 sec
    SL version: 7.0.1.66975
    created by: admin
       comment: Description for job Import PacBio DataSet
          tags:
          path: /global/seqfs/sdm/testing/pacbio/smrt_link5_1_denovo/smrtlink/userdata/jobs_root/003/003883

# Failed DMXing
...line 260, in version raise IOError("This does not appear to be a valid PacBio BAM file. Only datasets from RS II and Sequel instruments are supported by this program.")

# digging into BAM file
cd '/global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/subsample'
bam_n100='m64021_190510_094753_n100.subreads.bam'
sam_n100=$( echo $bam_n100| sed 's/bam/sam/')
$SMRT_BIN_701/bamsieve -n 100 \
    '/global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/PB1205_in/m64021_190510_094753.subreads.bam' \
    ./"$bam_n100"
samtools sort -n $bam_n100 |  samtools view -O SAM -h > $sam_n100

$SMRT_BIN_701/bamsieve -n 100 '/global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/PB1205_out_0625/PB1205_vanilla_004.subreads.bam' ./PB1205_vanilla_004.subreads.bam
#> IOError: This does not appear to be a valid PacBio BAM file. Only datasets from RS II and Sequel instruments are supported by this program

## repeat ICF on subsampled file
bam_in=$bam_n100
bam_out_white=$(echo $bam_n100 | sed 's/n100/n100_w/')
bam_out_black=$(echo $bam_n100 | sed 's/n100/n100_b/')
sam_n100_white=$( echo $bam_out_white| sed 's/bam/sam/')

date; time icecreamfinder.sh \
    in=$bam_in out=$bam_out_white outb=$bam_out_black \
    cq=f ow=t keepshortreads=f sambamba=f &
# completed
# Waiting on header to be read from a sam file.
# Time:                          1.548 seconds.
# Reads Processed:         278   0.18k reads/sec
# Bases Processed:       1931k   1.25m bases/sec
# Reads Out:               275
# Bases Out:             1890k
# Reads Filtered:            3   1.079%
# Total ZMWs:              132
# Bad ZMWs:                  2   1.515%
# Missing Adapter:           3   2.273%
# Hidden Adapter:            1   0.758%
# Avg Score Ratio:       0.787

samtools sort -n $bam_out_white | samtools view -O SAM -h > $sam_n100_white
# validating with bamsieve
$SMRT_BIN_701/pbindex $bam_out_white
$SMRT_BIN_701/bamsieve -n 10 $bam_out_white dev_null.bam


# that resolved issues with the header (at least it stays the same)
# now trying small-er bam for import/dmxing/hgap
thisDP='/global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763'
thisInDP="$thisDP/PB1205_in"
thisOutDP="$thisDP/PB1205_out_0626"
cd $thisOutDP

bam_in='m64021_190510_094753.subreads.bam'
bam_in_FP=$thisInDP/$bam_in
bam_in_sub15=$(echo $bam_in | sed 's/094753/001500/')
bam_in_sub15_FP=$thisOutDP/$bam_in_sub15
bam_out_white=$(echo $bam_in_sub15 | sed 's/000015/001501/')
bam_out_white_FP=$thisOutDP/$bam_out_white
bam_out_black=$(echo $bam_in_sub15 | sed 's/000015/001510/')
bam_out_black_FP=$thisOutDP/$bam_out_black

$SMRT_BIN_701/bamsieve --percent 15 \
    $bam_in_FP $bam_in_sub15_FP
#> completed    

date; time icecreamfinder.sh \
    in=$bam_in_sub15_FP out=$bam_out_white_FP outb=$bam_out_black_FP \
    cq=f ow=t keepshortreads=f sambamba=f jni=t &

# Found samtools 1.9
# Waiting on header to be read from a sam file.
# Time:                          736.265 seconds.
# Reads Processed:       1705k   2.32k reads/sec
# Bases Processed:      11610m   15.77m bases/sec
# Reads Out:             1670k
# Bases Out:            11268m
# Reads Filtered:        34979   2.051%
# Total ZMWs:           761035
# Bad ZMWs:               8381   1.101%
# Missing Adapter:        4107   0.540%
# Hidden Adapter:         4274   0.562%
# Avg Score Ratio:       0.758
# Score Cutoff:          0.640

# real 12m18.974s
# user 157m21.307s
# sys 4m8.518s


new_ds_name='PB1205_vanilla_005'
# Need 'original' xml >
$SMRT_BIN_701/dataset create \
  --type SubreadSet \
  --name $new_ds_name \
  --metadata "/global/seqfs/sdm/prod/pacbio/sequel/runs/64021/r64021_20190507_204150/4_A04/m64021_190510_094753.subreadset.xml" \
  --generateIndices \
    $new_ds_name.subreadset.xml \
    $bam_out_white_FP 

# validate 
$SMRT_BIN_701/dataset validate $new_ds_name.subreadset.xml
#> PB1205_vanilla_005.subreadset.xml is valid DataSet XML with valid ResourceId references

# import
SMRT_HOST_701='dint01.nersc.gov'  # 'https://dint01.nersc.gov:8243'
SMRT_PORT='9091'
XML2IMP=$thisOutDP/$new_ds_name.subreadset.xml
# CMD
$SMRT_BIN_701/pbservice import-dataset \
  --host $SMRT_HOST_701 --port $SMRT_PORT -u admin -p admin \
  "$XML2IMP"
  
# *DATASET SUMMARY*
          # id: 4674
        # uuid: e37bdd1d-9b92-594c-ba99-c807319c46d3
        # name: PB1205_vanilla_005
  # numRecords: 1670713
 # totalLength: 11268930663
   # createdAt: 2019-06-26T17:15:46.000Z
   # updatedAt: 2019-06-26T17:15:46.000Z
  # importedAt: 2019-06-27T00:21:18.577Z
   # createdBy: admin
        # tags:
    # isActive: true
 # parent UUID:
 # numchildren: 0
       # jobId: 3885
   # projectId: 1
        # path: /global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/PB1205_out_0626/PB1205_vanilla_005.subreadset.xml


# *JOB SUMMARY*
            # id: 3885
          # uuid: 6f438e8c-5940-4955-9bea-d7f0789bf9dd
          # name: import-dataset
         # state: SUCCESSFUL
    # project id: 1
     # jobTypeId: import-dataset
     # is active: true
     # createdAt: 2019-06-27T00:21:17.427Z
     # updatedAt: 2019-06-27T00:21:18.590Z
  # jobUpdatedAt: 2019-06-27T00:21:18.590Z
      # run time: 2 sec
    # SL version: 7.0.1.66975
    # created by: admin
       # comment: Description for job Import PacBio DataSet
          # tags:
          # path: /global/seqfs/sdm/testing/pacbio/smrt_link5_1_denovo/smrtlink/userdata/jobs_root/003/003885

# 2019-06-27
$SMRT_BIN_701/pbvalidate --type BAM --index  m64021_190510_001501.subreads.bam
# sh: /global/seqfs/sdm/testing/pacbio/smrt_link5_1_denovo/smrtlink/install/smrtlink-release_7.0.1.66975/bundles/smrttools/install/smrttools-release_7.0.1.66768/private/thirdparty/python/python_2.7.9/binwrap/../../../../thirdparty/readline/readline_6.3/lib/libreadline.so.6: no version information available (required by sh)

# - hard to debug issue with bash@nersc, according to Chris B
# tried: {cori,dint}, {bash,zsh}, {SMRT 6, SMRT 7} - no effect (failing)

# container
# Hi Alex, Here is an example of running hgap in shifter:
# $ shifter --image=bryce911/smrtlink:7.0.1.66768 \
# /smrtlink/bundles/smrttools/install/smrttools-release_7.0.1.66768/smrtcmds/bin/pbsmrtpipe \
# pipeline-id pbsmrtpipe.pipelines.polished_falcon_fat \
# --preset-json /global/projectb/scratch/brycef/pbmid/GNYAX/config/preset.json \
# --preset-json /global/projectb/scratch/brycef/pbmid/GNYAX/config/hgap_settings.json \
# -e eid_subread:/global/projectb/scratch/brycef/pbmid/GNYAX/config/input_GNYAX.xml \
# -o /global/projectb/scratch/brycef/pbmid/GNYAX/hgap

# To just get into the container you can use:
# $ shifter --image=bryce911/smrtlink:7.0.1.66768
# then to get to the root folder of the container:
# sh-4.2$ cd /

sshi
shifter --image=bryce911/smrtlink:7.0.1.66768 # stepping into interactive container mode
:sh-4.2$
SMRT_7C_BIN='/smrtlink/bundles/smrttools/install/smrttools-release_7.0.1.66768/smrtcmds/bin'
$SMRT_7C_BIN/pbvalidate # normal help, without 'libreadline.so.6: no version information available (required by sh)'
# native PacBio bam, subsampled with pbsieve
$SMRT_7C_BIN/pbvalidate --type BAM --index m64021_190510_001500.subreads.bam
# no error/warnings, so pre-ICF bam = validated
# ICF-ed PacBio bam
$SMRT_7C_BIN/pbvalidate --type BAM --index m64021_190510_001501.subreads.bam 
# no error/warnings, so post-ICF bam = validated
# now let's re-do dataset, within the container
new_ds_name='PB1205_vanilla_006'
$SMRT_7C_BIN/dataset create \
  --type SubreadSet \
  --name $new_ds_name \
  --metadata "/global/seqfs/sdm/prod/pacbio/sequel/runs/64021/r64021_20190507_204150/4_A04/m64021_190510_094753.subreadset.xml" \
  --generateIndices \
  $new_ds_name.subreadset.xml \
  'm64021_190510_001501.subreads.bam' 
# checking new version of dataset.xml (made without any err/warn :)
$SMRT_7C_BIN/pbvalidate --type SubreadSet PB1205_vanilla_006.subreadset.xml
# diff-ed: PB1205_out_0626/PB1205_vanilla_005.subreadset.xml : PB1205_out_0626/PB1205_vanilla_006.subreadset.xml
# only difference is time stamps, and UIDs = OK, no issues with how data were formatted

# import
SMRT_HOST_701='dint01.nersc.gov'  # 'https://dint01.nersc.gov:8243'
SMRT_PORT='9091'
XML2IMP=$(pwd)/$new_ds_name.subreadset.xml
$SMRT_7C_BIN/pbservice import-dataset \
  --host $SMRT_HOST_701 --port $SMRT_PORT -u admin -p admin \
  "$XML2IMP"
# success:
# *DATASET SUMMARY* 
#           id: 4691
#         uuid: 989035b4-4af4-ef69-426b-68c0495f2b3f
#         name: PB1205_vanilla_006
#   numRecords: 1670713
#  totalLength: 11268930663
#    createdAt: 2019-06-27T15:20:53.000Z
#    updatedAt: 2019-06-27T15:20:53.000Z
#   importedAt: 2019-06-27T22:37:49.366Z
#    createdBy: admin
#         tags: 
#     isActive: true
#  parent UUID: 
#  numchildren: 0
#        jobId: 3895
#    projectId: 1
#         path: /global/projectb/sandbox/rqc/ayspunde/work@sandbox/jira/seqqc-jira/seqqc-13763/PB1205_out_0626/PB1205_vanilla_006.subreadset.xml

# one more try: sorting ICFed bam
# samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] -T out.prefix [-@ threads] [in.bam]

samtools sort -o m64021_190510_001502.subreads.bam -O bam -n -T temp_sort -@ 16 m64021_190510_001501.subreads.bam

shifter --image=bryce911/smrtlink:7.0.1.66768 # stepping into interactive container mode
:sh-4.2$
export SMRT_7C_BIN='/smrtlink/bundles/smrttools/install/smrttools-release_7.0.1.66768/smrtcmds/bin'
    # ICF-ed PacBio bam
    # $SMRT_7C_BIN/pbvalidate --type BAM --index m64021_190510_001501.subreads.bam 
    # no error/warnings, so post-ICF bam = validated
    # now let's re-do dataset, within the container
new_ds_name='PB1205_vanilla_007'
$SMRT_7C_BIN/dataset create \
  --type SubreadSet \
  --name $new_ds_name \
  --metadata "/global/seqfs/sdm/prod/pacbio/sequel/runs/64021/r64021_20190507_204150/4_A04/m64021_190510_094753.subreadset.xml" \
  --generateIndices \
  $new_ds_name.subreadset.xml \
  'm64021_190510_001502.subreads.bam' &
      # checking new version of dataset.xml (made without any err/warn :)
      # $SMRT_7C_BIN/pbvalidate --type SubreadSet PB1205_vanilla_006.subreadset.xml
      # diff-ed: PB1205_out_0626/PB1205_vanilla_005.subreadset.xml : PB1205_out_0626/PB1205_vanilla_006.subreadset.xml
      # only difference is time stamps, and UIDs = OK, no issues with how data were formatted

# import
SMRT_HOST_701='dint01.nersc.gov'  # 'https://dint01.nersc.gov:8243'
SMRT_PORT='9091'
XML2IMP=$(pwd)/$new_ds_name.subreadset.xml
$SMRT_7C_BIN/pbservice import-dataset \
  --host $SMRT_HOST_701 --port $SMRT_PORT -u admin -p admin \
  "$XML2IMP"
# *JOB SUMMARY* 
#             id: 3898
#           uuid: fb2c1ccc-e064-4423-b56c-2f333f3eb204
#           name: import-dataset
#          state: SUCCESSFUL
