|
3 | 3 | import sys |
4 | 4 | import numpy as np |
5 | 5 |
|
| 6 | +# Parse the options from command line |
6 | 7 | try: |
7 | 8 | opts, args = getopt.getopt(sys.argv[1:],"i:n:k:t:",["idir=","njobs=","hacked-kallisto-path=","reference-transcriptome="]) |
8 | 9 | except getopt.GetoptError: |
9 | 10 | print ("getopterrror") |
10 | 11 | print ('usage is : \n python Trapnell_wrapper.py -i input_SRA_dir -k path-to-hacked-kallisto -t path-to-human-reference-transcriptome [-n number-of-processes-to-use]') |
11 | 12 | sys.exit(1) |
12 | | - |
| 13 | + |
| 14 | +# The four options used in this pipeline are |
13 | 15 | SRA_dir='' |
14 | 16 | num_proc=1 |
15 | 17 | kallipso_path='' |
|
25 | 27 | elif opt in ("-t","--reference-transcriptome"): |
26 | 28 | ref_transcriptome=arg |
27 | 29 |
|
28 | | - |
| 30 | +# Error message if function not called correctly |
29 | 31 | if (not SRA_dir) or (not kallipso_path) or (not ref_transcriptome): |
30 | 32 | print ('usage is : \n python Trapnell_wrapper.py -i input_SRA_dir -k path-to-hacked-kallisto -t path-to-human-reference-transcriptome [-n number-of-processes-to-use]') |
31 | 33 | sys.exit(1) |
32 | 34 |
|
| 35 | +# Execute process_SRA, which converts SRA to fastq |
33 | 36 | print('Extracting reads from SRAs...') |
34 | 37 | os.system('mkdir -p ./reads/') |
35 | 38 | os.system('rm -f ./reads/*') |
36 | 39 | os.system('python process_SRA.py -i '+SRA_dir+' -o ./reads/ -n '+str(num_proc)) |
37 | 40 |
|
| 41 | +# This removes files that Trapnell didn't use. (Specific for analysis of this dataset) |
38 | 42 | print('Removing files that Trapnell throws...') |
39 | 43 | relevant_fls=np.loadtxt('Files_to_keep.txt',dtype=str) |
40 | 44 | fls_to_remove=[cells for cells in os.listdir('./reads/') if cells.split('_')[0] not in relevant_fls ] |
41 | 45 | for fl in fls_to_remove: |
42 | 46 | os.system('rm '+'./reads/'+fl) |
43 | 47 |
|
44 | | - |
| 48 | +# Use kallisto to build index. should return same results as hacked kallisto |
45 | 49 | print('Generating the Kallisto index (with hacked kallisto)...') |
46 | 50 | os.system('mkdir -p ./kallisto_index') |
47 | 51 | os.system('rm -f ./kallisto_index/*') |
|
52 | 56 | num_ec = sum(1 for line in open('./kallisto_index/Trapnell_index.idx_ecmap.txt')) |
53 | 57 | print(num_ec) |
54 | 58 |
|
| 59 | +# Use hacked kallisto to generate TCCs (output is .class file that contains counts of ECs) |
55 | 60 | print('Generating TCC (with hacked kallisto)...') |
56 | 61 | os.system('mkdir -p ./transcript_compatibility_counts/') |
57 | 62 | os.system('rm -f ./transcript_compatibility_counts/*') |
58 | 63 | os.system('python get_pseudoalignments_paired_end.py -i ./reads/ -o ./transcript_compatibility_counts/ -k '+kallipso_path+ ' -t '+ index_path +' -n '+ str(num_proc)) |
59 | 64 |
|
| 65 | +# Process the .class file into a matrix. normalize such that each row sums to 1 |
60 | 66 | print('Generating TCC distribution...') |
61 | | - |
62 | 67 | os.system('python get_tcc_dist.py -i ./transcript_compatibility_counts/ -m '+str(num_ec)+' -t ./Trapnell_TCC.dat -d ./Trapnell_TCC_distribution.dat') |
63 | 68 |
|
| 69 | +# Process the TCC distribution matrix to get pairwise distributions |
64 | 70 | print('Generating pairwise distances...') |
65 | | - |
66 | 71 | os.system('python get_pairwise_distances.py ./Trapnell_TCC_distribution.dat ./Trapnell_TCC_pairwise_distance.dat '+str(num_proc)) |
0 commit comments