OSCR

Elevated intron retention implicates neuroinflammation in brains of individuals with alcohol use disorder.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § Materials and methods › Identification of IR events, quantification of intron content, and statistical analysis ↔ scripts/intron_calling.py, lines 62–142 · score 0.60 · HTSeq, uniquely mapped, interval, filtered, exons, events

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 496 lines · 25 KB · no license · 1 match

  1. import argparse;
  2. import os,shutil;
  3. import subprocess;
  4. import HTSeq,re,time;
  5. import collections;
  6. import pandas as pd;
  7. from multiprocessing import Pool;
  8. #
  9. parser = argparse.ArgumentParser()
  10. parser.add_argument('-b', metavar = 'input', dest='bamFile', help='Give input file fullname');
  11. parser.add_argument('-g', metavar = 'gtffile', dest='gtfFile', help='Give parsed gtf file fullname');
  12. parser.add_argument('-t', dest='thread', type=int, help='number of multiple thread number,>0');
  13. parser.add_argument('-o', metavar = 'output', dest='outputFile',help='Give output file fullname');
  14. parser.add_argument('-normal', metavar = 'common retain intron list', dest='norm_list',help='provide intron to be removed');
  15. parser.add_argument('-j', metavar = 'anchor length', type=int, dest='junction',help='Give a specific anchor length');
  16. parser.add_argument('-c', metavar = 'intronRead filter', dest='intron_read', type=int,help='Set a intron read filter');
  17. parser.add_argument('-p', metavar = 'intronPSI filter', dest='intron_psi', type=float,help='threshhold of psi value filter');
  18. parser.add_argument('-n', metavar = 'novel intron read read filter', dest='novel_read', type=int,help='Set a novel intron read filter threshhold');
  20. args = parser.parse_args();
  21. bamFile = args.bamFile;
  22. gtf = args.gtfFile;
  23. thread = args.thread;
  24. output = args.outputFile;
  25. norm_list = args.norm_list;
  26. junction = args.junction;
  27. intron_read_filter = args.intron_read;
  28. intron_psi_filter = args.intron_psi;
  29. novel_read_filter = args.novel_read;
  31. cpus= os.cpu_count();# get all available cpu counts
  32. if thread:
  33. thread=thread;
  34. thread= min(thread, cpus); # make sure the maximum thread not over available cpus!
  35. else:
  36. thread=1;
  38. if norm_list:
  39. normList=[line.strip() for line in open(norm_list, 'r')]; # read from provided remove list
  40. else:
  41. normList=[];
  43. if junction: # the anchor length setting
  44. junction= junction;
  45. else:
  46. junction= 10;
  48. if intron_read_filter: # the reads fall in intron region filter
  49. intron_read_filter = intron_read_filter;
  50. else:
  51. intron_read_filter = 10;
  53. if intron_psi_filter: # percentage of intron retained rate!
  54. intron_psi_filter = intron_psi_filter;
  55. else:
  56. intron_psi_filter = 0.05;
  58. if novel_read_filter: # the anchor length setting
  59. novel_read_filter = novel_read_filter;
  60. else:
  61. novel_read_filter= 15;
  62. #
  63. #
  64. def no_parallel_run():
  65. global features;
  66. global bamfile;
  68. counts = collections.Counter()
  69. bam = HTSeq.BAM_Reader(bamfile)
  70. for read in bam:
  71. #filter1, uniq mapp filter:the option as NM/NH/CC/CP/HI
  72. NH_list= [x for x in read.optional_fields if x[0]=='NH'];
  73. if len(NH_list)==0:
  74. NH_value= 10; # if not NH tag find assign as 10;
  75. elif len(NH_list)==1:
  76. NH_value = NH_list[0][1];
  77. if NH_value == 1 or read.aQual>=30: #filter1, uniq mapp filter
  78. cigar_list= [cstring.type for cstring in read.cigar]
  79. cigar_check= [e for e in cigar_list if e in ['D','I','S','H','P','X','=' ]];
  80. if len(cigar_check)==0:
  81. counts['uniq_mapped_read'] += 1;#count every unique mapped read for calculate rpkm value
  82. if cigar_list.count('N') ==0:# call directly
  83. gene_ids = []; iv_lens=[]
  84. for iv, val in features[ read.iv ].steps():
  85. if len(val) ==1:
  86. gene_ids.append(list(val)[0]);
  87. ivList= re.split('[:[,)/]',str(iv));
  88. iv_lens.append(str(int(ivList[3]) - int(ivList[2])))
  89. if len(gene_ids)>0 and len(set([x[:15] for x in gene_ids]))==1: # remove ambigous mapped read-map to two diff genes
  90. if len(gene_ids)==1:
  91. gene_id = gene_ids[0];
  92. counts[gene_id] += 1;
  93. elif len(gene_ids)>1: #strict the intron at least have length large than anchor
  94. for g in range(len(gene_ids)):
  95. gene_id = gene_ids[g];
  96. iv_len = int(iv_lens[g]);
  97. geneInfo = re.split('-|@', gene_id)[1]
  98. if(int(geneInfo)%2==1):
  99. counts[gene_id] += 1;
  100. elif (int(geneInfo)%2==0) and (iv_len > junction):
  101. counts[gene_id] += 1;
  102. elif cigar_list.count('N') ==1:
  103. gene_ids = []; iv_lens=[]
  104. for cigop in read.cigar:
  105. if cigop.type == "M":
  106. for iv, val in features[cigop.ref_iv].steps():
  107. if len(val) ==1:
  108. gene_ids.append(list(val)[0]);
  109. ivList= re.split('[:[,)/]',str(iv));
  110. iv_lens.append(str(int(ivList[3]) - int(ivList[2])))
  111. if len(set([x[:15] for x in gene_ids]))==1: # remove ambigous mapped read-map to two diff genes
  112. # get two terminal end of skip region
  113. if len(gene_ids)>1:
  114. skip_iv= re.split('[:[,)/]',str(read.cigar[1].ref_iv));
  115. new_query= HTSeq.GenomicInterval(skip_iv[0], int(skip_iv[2])-2,int(skip_iv[3])+3, skip_iv[5]); # make sure the enlarged map window to clarify the mapping region/s!
  116. query_ids=set()
  117. for iv, val in features[new_query].steps():
  118. query_ids |= val;
  119. if len(query_ids)>0:
  120. query_list= [int(re.split('-|@', x)[1]) for x in list(query_ids)];
  121. query_list.sort();# the mapped region list check
  122. if query_list[0]%2!=0 and query_list[-1]%2!=0:# if both end at exons
  123. for g in range(len(gene_ids)):
  124. gene_id = gene_ids[g];
  125. iv_len = int(iv_lens[g]);
  126. geneInfo = re.split('-|@', gene_id)[1]
  127. if(int(geneInfo)%2==1):
  128. counts[gene_id] += 1;
  129. elif (int(geneInfo)%2==0) and (iv_len > junction):
  130. counts[gene_id] += 1;
  131. else: # all should be add as novel intron/partial retained events
  132. novel_id= (list(query_ids)[0])[0:15] + '@' + skip_iv[0] + ':u' + skip_iv[2] + '-d' + skip_iv[3];
  133. counts[novel_id] += 1;
  134. elif len(gene_ids)==1: # if disruption occurs inner one feature
  135. skip_iv= re.split('[:[,)/]',str(read.cigar[1].ref_iv));
  136. gene_id = list(gene_ids)[0];
  137. geneInfo = re.split('-|@', gene_id)[1]
  138. if int(geneInfo)%2==0: # if read mapped within a intron, count it as novel events
  139. novel_id= gene_id[0:15] + '@' + skip_iv[0] + ':u' + skip_iv[2] + '-d' + skip_iv[3];
  140. counts[novel_id] += 1;
  141. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": HTSeq count finished!")
  142. return counts
  144. def do_parallel_run(part_N):
  146. global thread;
  147. global features;
  148. global tmp_prefix;
  150. bamfile_temp= tmp_prefix + "{0:0=2d}".format(part_N) + '.bam'
  151. counts = collections.Counter()
  152. bam = HTSeq.BAM_Reader(bamfile_temp)
  153. for read in bam:
  154. #filter1, uniq mapp filter:the option as NM/NH/CC/CP/HI
  155. NH_list= [x for x in read.optional_fields if x[0]=='NH'];
  156. if len(NH_list)==0:
  157. NH_value= 10; # if not NH tag find assign as 10;
  158. elif len(NH_list)==1:
  159. NH_value = NH_list[0][1];
  160. if NH_value == 1 or read.aQual>=30: #filter1, uniq mapp filter
  161. cigar_list= [cstring.type for cstring in read.cigar]
  162. cigar_check= [e for e in cigar_list if e in ['D','I','S','H','P','X','=' ]];
  163. if len(cigar_check)==0:
  164. counts['uniq_mapped_read'] += 1;#count every unique mapped read for calculate rpkm value
  165. if cigar_list.count('N') ==0:# call directly
  166. gene_ids = []; iv_lens=[]
  167. for iv, val in features[ read.iv ].steps():
  168. if len(val) ==1:
  169. gene_ids.append(list(val)[0]);
  170. ivList= re.split('[:[,)/]',str(iv));
  171. iv_lens.append(str(int(ivList[3]) - int(ivList[2])))
  172. if len(gene_ids)>0 and len(set([x[:15] for x in gene_ids]))==1: # remove ambigous mapped read-map to two diff genes
  173. if len(gene_ids)==1:
  174. gene_id = gene_ids[0];
  175. counts[gene_id] += 1;
  176. elif len(gene_ids)>1: #strict the intron at least have length large than anchor
  177. for g in range(len(gene_ids)):
  178. gene_id = gene_ids[g];
  179. iv_len = int(iv_lens[g]);
  180. geneInfo = re.split('-|@', gene_id)[1]
  181. if(int(geneInfo)%2==1):
  182. counts[gene_id] += 1;
  183. elif (int(geneInfo)%2==0) and (iv_len > junction):
  184. counts[gene_id] += 1;
  185. elif cigar_list.count('N') ==1:
  186. gene_ids = []; iv_lens=[]
  187. for cigop in read.cigar:
  188. if cigop.type == "M":
  189. for iv, val in features[cigop.ref_iv].steps():
  190. if len(val) ==1:
  191. gene_ids.append(list(val)[0]);
  192. ivList= re.split('[:[,)/]',str(iv));
  193. iv_lens.append(str(int(ivList[3]) - int(ivList[2])))
  194. if len(set([x[:15] for x in gene_ids]))==1: # remove ambigous mapped read-map to two diff genes
  195. # get two terminal end of skip region
  196. if len(gene_ids)>1:
  197. skip_iv= re.split('[:[,)/]',str(read.cigar[1].ref_iv));
  198. new_query= HTSeq.GenomicInterval(skip_iv[0], int(skip_iv[2])-2,int(skip_iv[3])+3, skip_iv[5]); # make sure the enlarged map window to clarify the mapping region/s!
  199. query_ids=set()
  200. for iv, val in features[new_query].steps():
  201. query_ids |= val;
  202. if len(query_ids)>0:
  203. query_list= [int(re.split('-|@', x)[1]) for x in list(query_ids)];
  204. query_list.sort();# the mapped region list check
  205. if query_list[0]%2!=0 and query_list[-1]%2!=0:# if both end at exons
  206. for g in range(len(gene_ids)):
  207. gene_id = gene_ids[g];
  208. iv_len = int(iv_lens[g]);
  209. geneInfo = re.split('-|@', gene_id)[1]
  210. if(int(geneInfo)%2==1):
  211. counts[gene_id] += 1;
  212. elif (int(geneInfo)%2==0) and (iv_len > junction):
  213. counts[gene_id] += 1;
  214. else: # all should be add as novel intron/partial retained events
  215. novel_id= (list(query_ids)[0])[0:15] + '@' + skip_iv[0] + ':u' + skip_iv[2] + '-d' + skip_iv[3];
  216. counts[novel_id] += 1;
  217. elif len(gene_ids)==1: # if disruption occurs inner one feature
  218. skip_iv= re.split('[:[,)/]',str(read.cigar[1].ref_iv));
  219. gene_id = list(gene_ids)[0];
  220. geneInfo = re.split('-|@', gene_id)[1]
  221. if int(geneInfo)%2==0: # if read mapped within a intron, count it as novel events
  222. novel_id= gene_id[0:15] + '@' + skip_iv[0] + ':u' + skip_iv[2] + '-d' + skip_iv[3];
  223. counts[novel_id] += 1;
  224. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": Process " + str(part_N) + " finished!")
  225. return counts
  227. def single_thread_RPKM(counter):
  228. geneSum = collections.Counter()
  229. countSets = collections.Counter()
  230. for k, v in counter.items():
  231. countSets[k] +=v;
  232. geneSum[k[:15]] +=v;
  234. del countSets["uniq_mapped_read"];
  235. keys= [key for key, value in countSets.items()];
  236. rpkm_list = [];
  237. for key in keys:
  238. value = countSets[key];
  239. list=[key, value] + [None]*10;
  240. extron_list= re.split('-|@|:|-', key);
  241. list[3] = geneSum[extron_list[0]];
  242. if len(extron_list)==5:
  243. extron_len = int(extron_list[4]) - int(extron_list[3]) + 1;
  244. list[9] = value*1000/extron_len;
  245. list[4] = round(value*(10**9)/(geneSum['uniq_mapped_rea'] * extron_len), 4);
  246. if (int(extron_list[1])%2==0):
  247. list[2] = 'intron';
  248. uExon = extron_list[0] + '-' + str(int(extron_list[1]) -1);
  249. dExon = extron_list[0] + '-' + str(int(extron_list[1]) + 1);
  250. upperExon = [[key, value] for key, value in countSets.items() if uExon in key];
  251. downExon = [[key, value] for key, value in countSets.items() if dExon in key]
  252. if len(upperExon) ==1:
  253. upperList = re.split('-|@|:|-', upperExon[0][0]);
  254. list[6] = upperExon[0][1]; list[10] = int(upperList[4])-int(upperList[3])+1;
  255. if len( downExon) == 1:
  256. downList = re.split('-|@|:|-', downExon[0][0]);
  257. list[7] = downExon[0][1]; list[11] = int(downList[4])-int(downList[3])+1;
  258. elif (int(extron_list[1])%2==1):
  259. list[2] = 'exon';
  260. elif len(extron_list)==4:
  261. list[2] = 'novel-intron';
  263. rpkm_list.append(list)
  264. return rpkm_list
  266. def do_RPKM_run(key):
  267. global geneSum;
  268. global countSets;
  270. value = countSets[key];
  271. list=[key, value] + [None]*10;
  272. extron_list= re.split('-|@|:|-', key);
  273. list[3] = geneSum[extron_list[0]];
  274. if len(extron_list)==5:
  275. extron_len = int(extron_list[4]) - int(extron_list[3]) + 1;
  276. list[9] = value*1000/extron_len;
  277. list[4] = round(value*(10**9)/(geneSum['uniq_mapped_rea'] * extron_len), 4);
  278. if (int(extron_list[1])%2==0):
  279. list[2] = 'intron';
  280. uExon = extron_list[0] + '-' + str(int(extron_list[1]) -1);
  281. dExon = extron_list[0] + '-' + str(int(extron_list[1]) + 1);
  282. upperExon = [[key, value] for key, value in countSets.items() if uExon in key];
  283. downExon = [[key, value] for key, value in countSets.items() if dExon in key]
  284. if len(upperExon) ==1:
  285. upperList = re.split('-|@|:|-', upperExon[0][0]);
  286. list[6] = upperExon[0][1]; list[10] = int(upperList[4])-int(upperList[3])+1;
  287. if len( downExon) == 1:
  288. downList = re.split('-|@|:|-', downExon[0][0]);
  289. list[7] = downExon[0][1]; list[11] = int(downList[4])-int(downList[3])+1;
  290. elif (int(extron_list[1])%2==1):
  291. list[2] = 'exon';
  292. elif len(extron_list)==4:
  293. list[2] = 'novel-intron';
  294. return list
  296. def single_thread_TPM(rpkm_list):
  297. listNum = len(rpkm_list);
  298. tpmList = [x[9] for x in rpkm_list];
  299. tpmSum = sum([x for x in tpmList if x != None]);
  300. tpm_list = [];
  301. for i in range(listNum):
  302. list_tpm = rpkm_list[i];
  303. if list_tpm[2] == 'intron':
  304. extron_list= re.split('-|@|:|-', list_tpm[0]);
  305. extron_len = int(extron_list[4]) - int(extron_list[3]) + 1;
  306. tpm_val = list_tpm[1]*(10**9)/(extron_len * tpmSum);
  307. list_tpm[5]= round(tpm_val, 4);
  308. if list_tpm[6] !=None:
  309. tpm_upper = list_tpm[6]*(10**9)/(list_tpm[10]* tpmSum);
  310. else:
  311. tpm_upper = 0;
  312. if list_tpm[7] !=None:
  313. tpm_down = list_tpm[7]*(10**9)/(list_tpm[11]* tpmSum);
  314. else:
  315. tpm_down = 0;
  316. tpm_flake=[tpm_upper, tpm_down]
  318. if sum(tpm_flake)!=0:
  319. list_tpm[8] = round(tpm_val*2/sum(tpm_flake), 3);
  320. else:
  321. list_tpm = list_tpm;
  322. tpm_list.append(list_tpm);
  323. return tpm_list
  325. def do_TPM_run(index):
  326. list_tpm = rpkm_list[index];
  327. if list_tpm[2] == 'intron':
  328. extron_list= re.split('-|@|:|-', list_tpm[0]);
  329. extron_len = int(extron_list[4]) - int(extron_list[3]) + 1;
  330. tpm_val = list_tpm[1]*(10**9)/(extron_len * tpmSum);
  332. list_tpm[5]= round(tpm_val, 4);
  334. if list_tpm[6] !=None:
  335. tpm_upper = list_tpm[6]*(10**9)/(list_tpm[10]* tpmSum);
  336. else:
  337. tpm_upper = 0;
  338. if list_tpm[7] !=None:
  339. tpm_down = list_tpm[7]*(10**9)/(list_tpm[11]* tpmSum);
  340. else:
  341. tpm_down = 0;
  342. tpm_flake=[tpm_upper, tpm_down]
  344. if sum(tpm_flake)!=0:
  345. list_tpm[8] = round(tpm_val*2/sum(tpm_flake), 3);
  346. else:
  347. list_tpm = list_tpm
  348. return list_tpm
  351. if __name__ == '__main__':
  352. global features;
  353. features = HTSeq.GenomicArrayOfSets( "auto", stranded=True )
  354. for line in open(gtf):
  355. fields = line.split( "\t" );
  356. chrom= fields[0]; start=int(fields[2]); end= int(fields[3])+1; strand=fields[4];
  357. name= fields[5].split('.')[0] + '-' + str(fields[12]) + '@' + fields[0] + ':' + str(fields[2]) + '-' + str(fields[3]);
  358. iv = HTSeq.GenomicInterval(chrom, start, end, strand)
  359. features[ iv ] += name;
  361. global bamfile;
  362. bamfile = bamFile;
  364. dir= '/'.join(bamfile.split('/')[0:-1]);
  365. s1= bamfile.split('/')[-1];
  366. sample= s1.replace('.bam', '');
  368. if output:
  369. if output[-1] =='/':
  370. output = output;
  371. else:
  372. output = output + '/'
  373. else:
  374. output = dir + '/'; # defualt located at your input bam file
  376. if thread==1:
  377. print("No parallel mode, thread =",thread)
  378. results= no_parallel_run();
  380. with open(output + sample + '_rawCount.txt', 'w') as f:
  381. for tag, count in results.items():
  382. f.write('{}\t{}\n'.format(tag, count))
  383. rpkm_list= single_thread_RPKM(results);
  384. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": RPKM expression calculated!")
  385. tpm_list = single_thread_TPM(rpkm_list);
  386. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": TPM and PSI calculated!")
  388. tpm_list = [x[:9] for x in tpm_list if x[2] !='exon']
  389. colname = [['id','count','feature','fullGeneCount','rpkm','tpm','uppExonCount','downExonCount','psi']];
  390. data= colname + tpm_list;
  391. data = pd.DataFrame(data[1:],columns=data[0]);
  392. data.to_csv(output + sample + '_intron_calling_Rawresult.txt', index=False, sep='\t');
  394. intron_list = [x for x in tpm_list if x[2] =='intron' and x[8]!=None and float(x[1]) >= intron_read_filter and float(x[8]) >= intron_psi_filter]
  395. # remove common retained intron list from normal samples
  396. intron_list = [x for x in intron_list if x[0] not in normList]
  397. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": Filter " + str(len(normList)) + " commonly retained introns in normal tissue!")
  398. intron_list = colname + intron_list;
  399. intron_data = pd.DataFrame(intron_list[1:], columns=intron_list[0]);
  400. intron_data.to_csv(output + sample + '_intron_candidates.txt', index=False, sep='\t');
  402. novel_list = [x for x in tpm_list if x[2] =='novel-intron' and x[1] >= novel_read_filter and x[1]/x[3] > 0.001]
  403. novel_list = colname + novel_list;
  404. novel_data = pd.DataFrame(novel_list[1:], columns=novel_list[0]);
  405. novel_data.to_csv(output + sample + '_novel-intron_candidates.txt', index=False, sep='\t');
  407. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": All intron-calling job done!")
  408. #
  409. else:
  410. print("thread =",thread)
  412. if not os.path.exists(dir + '/tmp_' + sample):
  413. os.makedirs(dir + '/tmp_' + sample)# create a folder for speeding calculation
  415. global tmp_prefix;
  416. tmp_prefix = dir + '/tmp_' + sample + '/tmp_'
  418. p0=subprocess.Popen("samtools view -@ " + str(thread) + " -c " + bamfile, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate()[0];
  419. readNum= int(p0.decode('utf-8').replace('\n', ''));# get bam total read number
  420. split_num= int(readNum/thread + 50);# set a split thread to split big bam to samll ones
  421. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": There are " + str(readNum) + " reads in this BAM file!" )
  423. p1=subprocess.Popen("samtools view -H " + bamfile +" > " + tmp_prefix + "header; samtools view -@ " + str(thread) +" " + bamfile + " | split - "+ tmp_prefix + " --numeric-suffixes -l " + str(split_num) + " --filter='cat " + tmp_prefix + "header - | samtools view -@ " + str(thread) + " -b - > $FILE.bam'", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate();
  424. # build index for temp bam files
  425. tmpBAM_list= os.listdir(dir + '/tmp_' + sample);
  426. tmpBAM_list= [dir + '/tmp_' + sample + '/' + x for x in tmpBAM_list if not 'header' in x]
  427. for tmpBAM in tmpBAM_list:
  428. p1=subprocess.Popen("samtools index -@ " + str(thread) + " " + tmpBAM, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate();
  429. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": Split BAM file into " + str(thread) + " small BAM files!")
  431. pool=Pool(processes=thread);
  432. process = pool.map(do_parallel_run, range(thread)); #multiple process
  434. result_list = [result for result in process];
  436. results=collections.Counter();
  437. for result in result_list:
  438. results += result;
  439. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": It may take 2~10 minutes to get the final results")
  441. with open(output + sample + '_rawCount.txt', 'w') as f:
  442. for tag, count in results.items():
  443. f.write('{}\t{}\n'.format(tag, count))
  444. shutil.rmtree(dir + '/tmp_' + sample) # remnove the temp folder
  446. global geneSum, countSets;
  447. geneSum = collections.Counter()
  448. countSets = collections.Counter()
  449. for k, v in results.items():
  450. countSets[k] +=v;
  451. geneSum[k[:15]] +=v;
  453. del countSets["uniq_mapped_read"];
  454. keys= [key for key, value in countSets.items()]
  456. pool=Pool(processes=thread);
  457. rpkm_process = pool.map(do_RPKM_run, keys); #multiple process
  458. rpkm_list = [list for list in rpkm_process];
  459. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": RPKM expression calculated!")
  461. # calculate tpm sum and psi value from flaking exons everate tpm
  462. listNum = len(rpkm_list);
  463. tpmList = [x[9] for x in rpkm_list];
  465. global rpkmList, tpmSum;
  466. tpmSum = sum([x for x in tpmList if x != None]);
  467. rpkmList = rpkm_list;
  468. pool=Pool(processes=thread);
  469. tpm_process = pool.map(do_TPM_run, range(listNum)); #multiple process
  470. tpm_list = [list_tpm for list_tpm in tpm_process];
  471. pool.close()
  472. pool.join()
  473. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": TPM and PSI calculated!")
  475. tpm_list = [x[:9] for x in tpm_list if x[2] !='exon']
  476. colname = [['id','count','feature','fullGeneCount','rpkm','tpm','uppExonCount','downExonCount','psi']];
  477. data= colname + tpm_list;
  478. data = pd.DataFrame(data[1:],columns=data[0]);
  479. data.to_csv(output + sample + '_intron_calling_Rawresult.txt', index=False, sep='\t');
  481. intron_list = [x for x in tpm_list if x[2] =='intron' and x[8]!=None and float(x[1]) >= intron_read_filter and float(x[8]) >= intron_psi_filter and float(x[8]) <0.5 and x[6]!=None and x[7]!=None and x[6] >=10 and x[7] >=10 and float(x[4]) >= 1 and float(x[4]) < 100];
  482. # remove common retained intron list from normal samples
  483. intron_list = [x for x in intron_list if x[0] not in normList]
  484. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": Filter " + str(len(normList)) + " commonly retained introns in normal tissue!")
  485. intron_list = colname + intron_list;
  486. intron_data = pd.DataFrame(intron_list[1:], columns=intron_list[0]);
  487. intron_data.to_csv(output + sample + '_intron_candidates.txt', index=False, sep='\t');
  489. novel_list = [x for x in tpm_list if x[2] =='novel-intron' and x[1] >= novel_read_filter]
  490. novel_list = colname + novel_list;
  491. novel_data = pd.DataFrame(novel_list[1:], columns=novel_list[0]);
  492. novel_data.to_csv(output + sample + '_novel-intron_candidates.txt', index=False, sep='\t');
  494. print(time.strftime("%Y-%m-%d %H:%M:%S") + ": All intron-calling job done!")
  495. #
  496. #

intron_calling.py at commit faaa43e, no license · at the source

Overview

Authors: Rudong Li1,2, Jill L Reiter1,2, Season K Wyatt-Johnson3,4, Chuanpeng Dong1,2, Caine S Smith5,6, Nick Green1,2, Hongyu Gao1,2, Sheketha R Hauser4,7, Manav Kapoor8,9, Julia Stevens5,6, R Dayne Mayfield10,11, Alison Goate8,9, Yue Wang2, Howard J Edenberg2,12, Richard L Bell7, Greg Trevor Sutherland5,6, Randy Brutkiewicz3,4, Yunlong Liu1,2
ORCID iDs: Rudong Li, Randy Brutkiewicz, Yunlong Liu
  1. Center for Computational Biology and Bioinformatics (CCBB), Indiana University (IU) School of Medicine, Indianapolis, IN 46202, USA
  2. Department of Medical and Molecular Genetics (MMGE), IU School of Medicine, Indianapolis, IN 46202, USA
  3. Department of Microbiology and Immunology, IU School of Medicine, Indianapolis, IN 46202, USA
  4. Stark Neurosciences Research Institute, IU School of Medicine, Indianapolis, IN 46202, USA
  5. New South Wales (NSW) Brain Tissue Research Centre, University of Sydney, Sydney, NSW 2006, Australia
  6. Charles Perkins Centre and School of Medical Sciences, University of Sydney, Sydney, NSW 2006, Australia
  7. Department of Psychiatry, IU School of Medicine, Indianapolis, IN 46202, USA
  8. Ronald M. Loeb Center for Alzheimer’s Disease, Department of Neuroscience, Icahn School of Medicine at Mount Sinai, New York, NY 10029, USA
  9. Department of Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai, New York, NY 10029, USA
  10. Waggoner Center for Alcohol and Addiction Research (WCAAR), University of Texas at Austin, Austin, TX 78712, USA
  11. Department of Neuroscience, University of Texas at Austin, Austin, TX 78712, USA
  12. Department of Biochemistry, Molecular Biology, and Pharmacology, IU School of Medicine, Indianapolis, IN 4620, USA
Institutions: Indiana University Indianapolis (United States); Indiana University School of Medicine; The University of Sydney (Australia); Icahn School of Medicine at Mount Sinai (United States); The University of Texas at Austin (United States)
Journal: Brain communications, volume 8, issue 5, article fcag264
Dates: received 5 November 2025; accepted 4 June 2026; published online 29 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1093/braincomms/fcag264 · PMID 42683190 · PMCID PMC13532580 · OpenAlex W7171617883
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), histology / microscopy (modality), human (organism), rat (organism), other condition (population), cellular / molecular (subfield)
Methods: Statistics, Connectivity, fMRI & imaging
Keywords: alcohol use disorder, double-stranded RNA, intron retention, neuroinflammation, alcohol-preferring P rats
Topic: RNA Research and Splicing (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: NIAAA NIH HHS (R28 AA012725)
Citations: not cited yet (Europe PMC); 76 references in the paper

Abstract

Intron retention, a form of alternative RNA splicing, can occur as part of normal gene regulation or result from disruption of the splicing machinery. Retained introns can potentially form double-stranded RNA, activating innate immune sensors and inflammation. This mechanism has been implicated in cancer but has not been studied in neuropsychiatric diseases like alcohol use disorder. We systematically analysed transcriptome-wide intron retention events in post-mortem brain tissue from 142 individuals (66 with alcohol use disorder and 76 controls), encompassing 320 region-specific samples from the superior frontal cortex, nucleus accumbens, central nucleus and basolateral amygdala. Analyses were adjusted for demographic, technical and biological covariates. Validation was performed in alcohol-preferring (P) rats using long-read sequencing. In complementary experiments, immunofluorescent staining was used to detect double-stranded RNA in rat brain tissue, while single-cell RNA-sequencing was performed to test activation of double-stranded RNA-sensing pathways in human brains. Brains from individuals with alcohol use disorder showed significantly higher total intron retention compared with controls, independent of age, with females showing greater increases than males. A total of 368 introns were positively associated with alcohol use disorder, and these introns were significantly longer and had weaker splice acceptor sites compared with non-associated introns. Genes harbouring these intron retention events were enriched in Purkinje neurons, visual cortex neurons and oligodendrocytes. Computational predictions indicated these long introns could form duplex RNA structures. Increased double-stranded RNA was confirmed experimentally in multiple brain regions of alcohol-consuming rats, where it co-localized primarily with neuronal nuclei and dendrites. In individuals with alcohol use disorder, we found that multiple pathways including double-stranded RNA responses, neuroinflammation, interferon and NF-κB signalling, adaptive immunity and apoptosis were activated. In addition, NeuN-positive neuronal counts significantly decreased in both the prefrontal and visual cortices. Furthermore, single-cell analysis demonstrated upregulation of TICAM1, the target of double-stranded RNA sensor TLR3, in oligodendrocytes, as well as widespread activation of downstream inflammatory pathways across glial and neuronal cell types. These findings provide the first evidence that chronic alcohol consumption promotes an overall increase of intron retention in the brain and is associated with the presence of double-stranded RNA. Furthermore, the double-stranded RNA may contribute to neuronal loss and brain pathology by activating a neuroinflammatory response.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

cpdong/IntronNeoantigen

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: faaa43e2549522e0f225d2429885fc400e0e20d8, 2 August 2020
Languages: Python (4)
Size: 11 files, 4 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (4 files), NumPy (3 files), Biopython (2 files), BEDTools (1 file), SAMtools (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
5 files
At the source: github.com/cpdong/IntronNeoantigen

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 4 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

Human brain sample data are available via the NCBI BioProject database: BLA (PRJNA551909): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA551909, CE (PRJNA551908): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA551908, NAC (PRJNA551775): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA551775, SFC (PRJNA530758): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA530758 and PFC (PRJNA781630): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA781630. The Python tool IR-NeoAg is available via GitHub, https://github.com/cpdong/IntronNeoantigen. WebCSEA is available at https://bioinfo.uth.edu/webcsea/. The R package clusterProfiler is available via Bioconductor, https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html. SPOT-RNA is available via GitHub, https://github.com/jaswindersingh2/SPOT-RNA. The RNA secondary structure tool bpRNA is available via GitHub, https://github.com/hendrixlab/bpRNA. This study did not generate any new code.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 18 authors, 5 keywords, 1 funder, 70 references.

Cite

This paper

Li, R., Reiter, J. L., Wyatt-Johnson, S. K., Dong, C., Smith, C. S., Green, N., Gao, H., Hauser, S. R., Kapoor, M., Stevens, J., Mayfield, R. D., Goate, A., Wang, Y., Edenberg, H. J., Bell, R. L., Sutherland, G. T., Brutkiewicz, R., & Liu, Y. (2026). Elevated intron retention implicates neuroinflammation in brains of individuals with alcohol use disorder. Brain communications, 8(5), fcag264. https://doi.org/10.1093/braincomms/fcag264

BibTeX

@article{li2026elevated,
author = {Li, Rudong and Reiter, Jill L and Wyatt-Johnson, Season K and Dong, Chuanpeng and Smith, Caine S and Green, Nick and Gao, Hongyu and Hauser, Sheketha R and Kapoor, Manav and Stevens, Julia and Mayfield, R Dayne and Goate, Alison and Wang, Yue and Edenberg, Howard J and Bell, Richard L and Sutherland, Greg Trevor and Brutkiewicz, Randy and Liu, Yunlong},
title = {{Elevated intron retention implicates neuroinflammation in brains of individuals with alcohol use disorder}},
journal = {Brain communications},
year = {2026},
month = jul,
volume = {8},
number = {5},
pages = {fcag264},
publisher = {Oxford University Press},
issn = {2632-1297},
doi = {10.1093/braincomms/fcag264},
url = {https://doi.org/10.1093/braincomms/fcag264},
pmid = {42683190},
pmcid = {PMC13532580}
}

RIS

TY - JOUR
AU - Li, Rudong
AU - Reiter, Jill L
AU - Wyatt-Johnson, Season K
AU - Dong, Chuanpeng
AU - Smith, Caine S
AU - Green, Nick
AU - Gao, Hongyu
AU - Hauser, Sheketha R
AU - Kapoor, Manav
AU - Stevens, Julia
AU - Mayfield, R Dayne
AU - Goate, Alison
AU - Wang, Yue
AU - Edenberg, Howard J
AU - Bell, Richard L
AU - Sutherland, Greg Trevor
AU - Brutkiewicz, Randy
AU - Liu, Yunlong
TI - Elevated intron retention implicates neuroinflammation in brains of individuals with alcohol use disorder
T2 - Brain communications
J2 - Brain Commun
PY - 2026
DA - 2026/07/29
VL - 8
IS - 5
SP - fcag264
SN - 2632-1297
PB - Oxford University Press
DO - 10.1093/braincomms/fcag264
UR - https://doi.org/10.1093/braincomms/fcag264
LA - en
ER -

CSL-JSON

{
"id": "10.1093/braincomms/fcag264",
"type": "article-journal",
"title": "Elevated intron retention implicates neuroinflammation in brains of individuals with alcohol use disorder",
"container-title": "Brain communications",
"author": [
{
"family": "Li",
"given": "Rudong"
},
{
"family": "Reiter",
"given": "Jill L"
},
{
"family": "Wyatt-Johnson",
"given": "Season K"
},
{
"family": "Dong",
"given": "Chuanpeng"
},
{
"family": "Smith",
"given": "Caine S"
},
{
"family": "Green",
"given": "Nick"
},
{
"family": "Gao",
"given": "Hongyu"
},
{
"family": "Hauser",
"given": "Sheketha R"
},
{
"family": "Kapoor",
"given": "Manav"
},
{
"family": "Stevens",
"given": "Julia"
},
{
"family": "Mayfield",
"given": "R Dayne"
},
{
"family": "Goate",
"given": "Alison"
},
{
"family": "Wang",
"given": "Yue"
},
{
"family": "Edenberg",
"given": "Howard J"
},
{
"family": "Bell",
"given": "Richard L"
},
{
"family": "Sutherland",
"given": "Greg Trevor"
},
{
"family": "Brutkiewicz",
"given": "Randy"
},
{
"family": "Liu",
"given": "Yunlong"
}
],
"container-title-short": "Brain Commun",
"volume": "8",
"issue": "5",
"page": "fcag264",
"DOI": "10.1093/braincomms/fcag264",
"PMID": "42683190",
"PMCID": "PMC13532580",
"ISSN": "2632-1297",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/braincomms/fcag264",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
29
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: Biopython, BEDTools, SAMtools, 2 other tools, cellular / molecular, 2 references
[2] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: Biopython, BEDTools, SAMtools, 3 other tools, genetics / omics, 1 reference
[3] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: Biopython, BEDTools, SAMtools, 2 other tools, cellular / molecular, 2 references
[4] doi:10.1523/eneuro.0484-25.2026
Epigenetic and Transcriptomic Impacts of Ethanol Vary by Brain Region and Extent of Exposure.
Journal: eNeuro
In common: genetics / omics, 6 references
[5] doi:10.1038/s41592-026-03211-w [code]
Spatial isoform sequencing at single-cell resolution reveals cell-type-specific spatial isoform variability in multiple brain cell types.
Journal: Nature methods
In common: Biopython, BEDTools, SAMtools, 3 other tools, genetics / omics
[6] doi:10.1002/alz.71271 [code]
Intracellular protein GBF1 displays significant associations with amyloid pathology in Alzheimer's disease.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: Biopython, BEDTools, SAMtools, 2 other tools, histology / microscopy, cellular / molecular
[7] doi:10.1101/gr.280394.124 [code]
De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes.
Journal: Genome research
In common: Biopython, BEDTools, scikit-learn, 2 other tools, cellular / molecular, 2 references
[8] doi:10.1007/s00262-026-04390-3 [code]
Identification and prioritisation of tumour antigen candidates from 79 glioblastoma transcriptomes.
Journal: Cancer immunology, immunotherapy : CII
In common: BEDTools, SAMtools, scikit-learn, 2 other tools, genetics / omics, other condition, cellular / molecular, 1 reference
[9] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: BEDTools, SAMtools, scikit-learn, 2 other tools, genetics / omics, cellular / molecular, 1 reference
[10] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: BEDTools, SAMtools, scikit-learn, 2 other tools, genetics / omics, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.