bam file thred score 및 flag확인용

#!/bin/bash
#$ -S /bin/bash
#$ -N g1
#$ -cwd

for x in /scratch/sjoh/GBM/cancer/WXS_T/*.bam
do
    echo "$x"
    y=${x##*/}
    echo "$y"
    samtools view $x | cut -f2,11 > /scratch/sjoh/GBM/cancer/WXS_T/quality/$y.result
    echo "/scratch/sjoh/GBM/cancer/WXS_T/quality/$y.result exported"
done
rm tmp

#------------------
# This script is to investigate thred score in bam file
#-----------------
'''
Before running python script, I parsed bam file and extracted flag and thred score filed.
    flag : Information that a read is properly mapped, paired-end, and etc
    thred score : base quality\
#------------ script -------------------------
    #!/bin/bash
    #$ -S /bin/bash
    #$ -N g1
    #$ -cwd

    for x in /scratch/sjoh/GBM/cancer/WXS_T/*.bam
    do
        echo "$x"
        y=${x##*/}
        echo "$y"
        samtools view $x | cut -f2,11 > /scratch/sjoh/GBM/cancer/WXS_T/quality/$y.result
        echo "/scratch/sjoh/GBM/cancer/WXS_T/quality/$y.result exported"
    done
    rm tmp
#----------------------------------------------
linux cluster

module load Python_3.4.1

'''

from collections import Counter
import os, fnmatch, string, subprocess
import json
# List of files *.result which contain flag and thred score
fls=fnmatch.filter(os.listdir('.'),'*bam.result')

# Do process
for i in range(0,len(fls)):
    print(fls[i])
# thread score parsing
    #commands='cut -f2 '+fls[i]+' > tm2'
    commands='cut -f2 '+fls[i]
    #os.system(commands) # Invoke process to extract thread score part
    output=subprocess.check_output(commands,shell=True)
    tm3=output.decode('utf-8')
       
# Load the file and Get thred counts
    #f=open('tm2','r')
    #tm3=f.read()
    tm4=''.join(tm3.split('\n'))
    tm5=json.dumps(Counter(tm4))
    #f.close()
# Export the result
    new_fname=fls[i]+'.thred_counts'
    with open(new_fname,'w') as file:
        file.write(tm5)
       
# thread score parsing
    #commands='cut -f1 '+fls[i]+' > tm2'
    commands='cut -f1 '+fls[i]
    #os.system(commands) # Invoke process to extract flag part
    output=subprocess.check_output(commands,shell=True)
    tm3=output.decode('utf-8')
# Load the file and Get thred counts
    #f=open('tm2','r')
    #tm3=f.read()
    #tm3=output
    tm4=tm3.split('\n')
    tm5=json.dumps(Counter(tm4[0:len(tm4)-1]))
    #f.close()
# Export the result
    new_fname=fls[i]+'.flag_counts'
    with open(new_fname,'w') as file:
        file.write(tm5)