xinglab/pairadise

Paired Replicate Analysis of Allelic Differential Splicing Events

pairadise's Issues

shinyapp not avaliable

Hello, I cannot open the shinyapp website you provided in your passage. Is there any solution?

step 1 error: [qsub: illegal -c value]

sorry to bother you again,
when I performed step 1:

qsub -cwd -l h_vmem=90G qsub_files/pairadise_step1_sample1.qsub

resulted in:
[qsub: illegal -c value]

usage: qsub [-a date_time] [-A account_string] [-b secs]
[-c [ none | { enabled | periodic | shutdown |
depth= | dir= | interval=}... ]
[-C directive_prefix] [-d path] [-D path]
[-e path] [-h] [-I] [-j oe|eo|n] [-k {oe}] [-K ]
[-l resource_list] [-m n|{abe}] [-M user_list] [-N jobname] [-o path]
[-p priority] [-P proxy_user [-J <jobid]] [-q queue] [-r y|n]
[-S path] [-t number_to_submit] [-T type] [-u user_list] [-w] path
[-W additional_attributes] [-v variable_list] [-V ] [-x] [-X] [-z] [script]

I know less about qsub, what can I do to continue this process?

Thank you.

Why number of cores limited to 8 ?

Hi,

I am just curious why there is a limit of 8 cores in foreach function. I have access to cluster that have nodes with 40 cores, I am wondering if I could use more cores than 8.

PAIRADISE/pairadise/src/pairadise_model/R/pairadise.R

Lines 83 to 84 in 105b168

	if(numCluster <= 0 | numCluster > 8){
	stop("Error: Number of clusters must be between 1 and 8")

Thanks,
Goutham A

identifying the functional SNP

Hi,
I have performed the PAIRADISE analysis in our own data, and we detected some significant SNP-splicing event pairs, there are two questions that confused me:

1. Can I directly treated these SNPs as fuctional SNP?
2. The number of the pairs were limited, I guess it may be attributed to restricting the candidate(heterozygous) SNPs to the flanking constitutive exons region. As far as I know, the SNPs within the alternative exon can also affect the alternative splicing, is there any solution to include these SNPs as well?

Thank you.

some bugs

Hi,

First question, the manual states that 'RNA editing information' is an optional file, but many parts of the script use rnaedit.file. So, I vim pairadise_assign/scripts/assign.py:

fh.write("pairadise_assign -o " + project_name + "/" + fn + " -v genotype/" + fid + "/ --gz\n")

`def read_in_rna_editing(self):

sys.stdout.write('processing: ' + self.editFile + '\n')
e = defaultdict(lambda:defaultdict(list))
if self.editFile == '': return e  ##funong
fin = open(self.editFile)
header = []
firstline = True
for line in fin:
  if firstline:
    header = line.rstrip().split()
    firstline = False
  else:
    fields = line.rstrip().split()
    result = {}
    for i,col in enumerate(header):
      result[col] = fields[i]
    e[result['chromosome'][3:]][(int(result['position'])-1)/100].append(int(result['position'])-1)  

`

Second, in pairadise_assign/scripts/assign.py, I found some len(reads1[qname]) != len(reads2[qname]), so this step

if all([True if reads1[qname][i].pos==reads2[qname][i].pos else False for i in range(len(reads1[qname]))]): # reads have same starting position in hap1 and hap2

will report an index error. So, I vim pairadise_assign/scripts/assign.py:

`for qname in snpreads1[0]:

if qname in reads2:
    print("reads1 is:", reads1[qname], ",reads2[qname] is:", reads2[qname])##funong
    print("reads1 length is:", len(reads1[qname]), ",reads2 length is:", len(reads2[qname]))##funong
    if len(reads1[qname]) == len(reads2[qname]):##funong
        if all([True if reads1[qname][i].pos == reads2[qname][i].pos else False for i in range(len(reads1[qname]))]):  # reads have same starting position in hap1 and hap2
            if self.num_mismatches(reads1[qname][0]) < self.num_mismatches(reads2[qname][0]):
                spec1.append(qname)
    else:##funong
        print("wrong1:reads1 length is not equal reads2", ",reads1 length is:", len(reads1[qname]),",reads2 length is:", len(reads2[qname]))##funong
else:
    hap1only.append(qname)

`

After making the above modifications, it seems that the script can run, but I am not sure if it will have a significant impact on the results. Hope to receive your feedback.

personalize genome

Hi,
I wanted to know how can I change the -r and -e arguments before running qsub step 1. I see that after running run_preprocess.py several qsub files are generated for each sample. I want to change the -r genome/hg19.fa and -e genome/Human_AG_all_hg19_v2.txt in the qsub1 because I have my sequences aligned in hg38. I see there is this personalize.py script; however it is not clear for me how to use it.
Best,
Bea

personalize report and assignment error

Hello, I have two problem.

The first one, I ran the personalize procedure as:
pairadise_personalize -o output/sample name -v genotype/sample name -r genome/hg19.fa -e genome/Human_AG_all_hg19_v2.txt --rnaedit --gz, but the report.personalize.txt shows an outcome of # number of hap1 SNPs: 0, # number of editing sites: 2576459. Is that normal?

The second one is in the assignment procedure:
Traceback (most recent call last):
File ".../bin/pairadise_assign", line 8, in
sys.exit(main())
File ".../pairadise_assign/scripts/assign.py", line 451, in main
p.haplotype_assignment()
File ".../pairadise_assign/scripts/assign.py", line 332, in haplotype_assignment
os.rename(input_files1[0],os.path.join(self.outDir,"hap1.bam"))
IndexError: list index out of range
How can I solve this?

Thank you.

pairadise computes p-value problems

Hi Eric Kutschera,

When I extract the result file RI.MATS.JC.txt then extract the key information in it, as shown below:

RI.p<-RI[,c("ID", "IJC_SAMPLE_1", "SJC_SAMPLE_1", "IJC_SAMPLE_2",

    "SJC_SAMPLE_2", "IncFormLen", "SkipFormLen")]

But when I use pairadise to calculate the p-value. He had a problem with data mismatch.As shown below:

Error in load.data(my.data) :
Error with data: some data are not matched pairs.

I also looked at your scripts on github and found that the two sample sizes don't match up. But we're using rMATs to start with the bam file and he's able to compute p-values.

But I also checked the commands in your paired_model.R file and found that it was the same as the commands we ran. Why did yours have no errors and corresponding results while mine did not? Even though I filled 0 for it to make their sample numbers consistent, the final result was still different from that of rMATs.

So I take the liberty of writing to ask you. Sorry about that.

Here I first say thank you!

Could you please provide an example of the RNA edit file?

Due to some reason, the file containing RNA edit information cannot be downloaded from RADAR. Could you please tell me what is the data form of this txt, which will make it available for us to use the RNA edit file downloaded from other websites, such as the Rediportal. Or it would be appreciated if you put this RADAR file into github.

Thank you.