xxxxx

python调用系统命令执行#!/usr/bin/env python# -*- coding: utf-8 -*-"""这个脚本用来循环调取shell命令,同时读取程序的输入文件,输出内容到文件"""import subprocessreport = open('hisat2_report2.txt', 'w')with open('hisat2_list2.txt', 'r') as f

Extract User Defined Region From An Chromosome Fasta File小麦基因组一条染色体序列至少有500M，我们想要获得染色体某一段的序列还是比较难办的。除了samtool可以做到，今天发现一个用python写的脚本pyfaidx，使用起来也比较方便。首先安装pip install --user pyfaidx #安装在自己账户的 $HOME/.loc

from Bio import SeqIOname = open('1.txt', 'r')o = open('1.fa', 'r')wanted = set(line.rstrip("\n") for line in name)records = (r for r in SeqIO.parse('unigene_seq_5.fasta', "fasta") if r.id not in

RNA_Seq差异表达分析流程1、数据下载ftp.sra.ebi.ac.uk/vol1/fastq/SRR122/005/SRR1228245/SRR1228245_1.fastq.gz;ftp.sra.ebi.ac.uk/vol1/fastq/SRR122/005/SRR1228245/SRR1228245_2.fastq.gzftp.sra.ebi.ac.uk/vol1/fastq/SRR12

GFF3 TO GTFgff3格式是使用gmap软件得到的。 输入文件gff3的格式如下：chr1A IWGSCv1.0_gmap gene 11740 12074 . + . ID=TRIAE_CS42_1AS_TGACv1_023354_AA0082670.1.path1;Name=TRIAE_CS42_1AS_TGACv1_023354_AA0082670.1

fasta转成txt有的时候需要将fasta格式的序列文件转换成以tab键分割的txt文件写了一个简单的命令。用法如下：fasta2txt -i input.fa -o out.txt代码如下#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ = "Sheng-Wei Ma"import click #需要你的电脑安装click包@

计算核酸序列的GC含量 （GC content）#!/usr/bin/env python# -*- coding: utf-8 -*- __author__ = 'shengwei ma'__author_email__ = 'shengweima@icloud.com'from Bio import SeqIOfrom Bio.SeqUtils import GCfor rec in S

给GFF3格式文件添加fasta格式　　是不是没见过带有序列的gff3格式。为啥这么做，这就要说到我最近在做的东西了。Jbrowse是一款基因组可视化浏览器。可以将基因组可视化以及大部分以基因组为基础的可视化，比如reads、SNP、QTL、GWAS、gene。支持fasta,bam,vcf,gff3等格式文件。说了这么多，给个实例，自己慢慢体会。同时附上官网地址和Genome Biology上的论

看到一个现金求助的题目：http://www.timedoo.com/task-id-1194.html代码如下：#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ = "Sheng-Wei Ma"from Bio import SeqIOfrom collections import OrderedDictrecords = (r f

删除 setup.py 安装的 Python 软件包　　在 CentOS 5.5 上通过 setup.py 安装了一个软件包。删除的时候发现 setup.py 没有 uninstall 选项。增加 –record 参数重新安装软件包，执行命令： python ./setup.py install --record install.txt删除安装文件，执行命令： cat install.txt |

输入文件input_v1.0.txt （三列，分别是 *.1.fastq.gz,*2.fastq.gz , *.sam）hisat2运行参数与流程（hisat2_IWGSCv1.0.py）#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ = 'shengwei ma'__author_email__ = 'shengweima@icl

python 处理xml文件　　最近基因注释需要查阅文献是否报道过。由于基因很多，想了一个办法。NCBI上每个蛋白有关的登录号下会有文献的题目。根据序列比对结果，然后调取对应的文献。首先获取小麦族（147389）所有的199754条蛋白序列，截止日期是17-5-22.下载的格式是INSDSeq XML格式。下载之后需要转换成表格形式首先需要编辑下下载的xml文件，分别在文件头以及文件尾分别添加如下内

Mikado - pick your transcript: a pipeline to determine and select the best RNA-Seq predictionMikado is a lightweight Python3 pipeline to identify the most useful or “best” set of transcripts from multi

经常遇到提取问题，一旦要提取的文件很大，关键字很多，可以使用集合#!/usr/bin/env python# -*- coding: utf-8 -*- __author__ = 'Shengwei Ma'__author_email__ = 'shengweima@icloud.com'with open('3.txt', 'r') as f1: a = set(line.stri

基因表达量FPKM计算之后，可以根据其值分组。   输入格式：Geneid UN003324 UN004629 UN011750 UN013498 UN016818 UN018363 UN018537 UN019619 UN019632 UN026711 UN054416 UN056756 UN058273root_Z10

NCBI   EST 文库格式转换#!/usr/bin/env python# -*- coding: utf-8 -*-with open('1.txt', 'r') as f: a = [] b = [] for num, line in enumerate(f): if 'Lib' not in line: line1

1.txtUN227692 chr6B 558820383 558820604 . 100.0 99.1 216 2 0 4UN113387 chr7A 683635472 683635624 . 100.0 100.0 153 0 0 0UN128584 chr7D 27592786 27593326 . 100.0 100.0 541 0 0 0UN170802 chr4B 5053

##gff-version 3# Generated by GMAP version 2016-06-09 using call: gmapl.sse42 -D /export/data/ -d NRGenome --trim-end-exons=10 -t 32 --canonical-mode=2 --allow-close-indels=2 -B 4 -f 4 -n 0 ./unig

Sequence CleanerDescriptionI want to share my script using Biopython to clean sequences up. You should know that analyzing poor data takes CPU time and interpreting the results from poor data

负号表示反向序列，并且start and end 的位置是序列反向互补之后的位置Name strand start end total_lengthUN000226 -2 1088 1411 1431UN073473 3 861 1082 1158UN082299 -1 838 1014 1105UN064320 2 791 1006 1150UN070736 -3 780 1271

软件exonerate输出的结果如下，想要获得比对上的target序列Command line: [./exonerate INPUT/UN029382.fa INPUT/scaffold125532.fa --model est2genome --showtargetgff TRUE --showvulgar no --showalignment yes --alignmentwidth 2

>maker-lcl|TGACv1_scaffold_433174_5DL_UN109209-exonerate_est2genome-gene-0.0-mRNA-1 transcript offset:19 AED:0.00 eAED:0.00 QI:19|1|1|1|0|0|3|270|77GGATGTTGCTACTTGCTAGATGGAAATGGAAACATGCAAGTTGGAACTTCC

根据blast结果选择完全overlap的序列，此时没有考虑identity的高低Query Target overlap_length identity Query_length Target_length Percent Q_start Q_end T_start T_end strandUN065663 maker-lcl|TGACv1_scaffold_643943_U_UN14

chr1AL-3889054 2870 3970 + Traes_1AL_0007E2761.1chr1AL-3889054 3220 3970 + Traes_1AL_0007E2761.2chr1AL-3877729 317 552 + Traes_1AL_002FAE6E8.1chr1AL-3947740 1 401 - Traes_1AL_00485097A.1chr1AL-394

cat 1.txta1 = (DESCRIPTION = (ADDRESS = (PROTOCOL = TCP)(HOST = oracle)(PORT = 1000)) (CONNECT_DATA = (SERVER = aaa) (SERVICE_NAME = aaa)a2 = (DESCRIPTION = (ADD

文本内容如下：192.168.100.125 UNKNOWN w0100441 [03/Jun/2015:16:16:26 +0800] 226 119096 0.014..........................63.54.78.89 UNKNOWN w0500465 [03/Sep/2015:16:16:26 +0800] 275 119098 0.01

根据GFF3文件统计外显子大小和数量以及内含子大小#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ = "Sheng-Wei Ma"with open('TGACv1.cdna.gff3', 'r') as f: for line in f: lin = line.strip().split('\t')