xxxxx

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包@

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

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

两个文件的合并分别有文件一和文件二，文件一和文件二都有一列内容是相同的。根据此列内容可将两者结合起来。文件一的格式如下 AEGTA04333 LOC_Os01g33520|LOC_Os01g33520.1 2.70E-25 LOC_Os01g33520|LOC_Os01g33520.1 PF02902 Ulp1 protease family, C-terminal cat

推荐：蛋白在线比对网站https://www.ebi.ac.uk/Tools/hmmer/提交一条蛋白序列去调取数据库中库相似的蛋白序列。结果与blastp结果类似。但是比blastp结果多出来很多其他的东西，比如，结构域检测，序列的多重比对.结果部分共分为四项，Score, Taxonomy, Domain,Download.

​​​### 遗传图谱比较或者遗传图和物理图谱比较今天在Bioinformatics杂志上看到一篇关于遗传图谱之间比较的文章，感觉不错，特介绍如下。题目是：The Genetic Map Comparator: a user-friendly application to display and compare genetic mapsPMID: 28073759 DOI: 10.10

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

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

绘制基因位于染色体上的图（map gene to chromosome）最近有一批基因要根据其在染色体上的位置将其展示出来，想到的一种解决办法是：参照遗传连锁图的画法将其展示出来，所以能绘制遗传连锁图的软件都可以使用，如mapchart，genetic-mapper等。另外一种方法就是寻找合适的绘制软件，运气不错，还真让我找到一个在线网站可以绘制此图（http://mg2c.iask.in/mg2

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

PortcullisPortcullis stands for PORTable CULLing of Invalid Splice junctions from pre-aligned RNA-seq data. It is known that RNAseq mapping tools generate many invalid junction predictions, particularl

Carthagene安装 最近在写文章，需要用carthagene做遗传连锁图。记录下安装过程。 下载地址： http://mulcyber.toulouse.inra.fr/frs/download.php/1199/carthagene-1.3.beta-Linux-x86_64.sh 我的是centos6，所以我下载的是x86_64 当然安装之前需要安装一些依赖包: • tcl

计算核酸序列的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

detectMITE 安装以及使用　　安装之前需要安装两个其他软件。第一个就是matlab,第二个是cd-hist.1、matlab的安装下载好之后：第一步是挂载iso文件mount -o loop Matlab801_MacUnix.iso /mnt/第二步是安装cd /mnt/./install -mode silent -agreeToLicense yes -fileInstallatio

makeblastdb注意事项本地blast，需要自行构建数据库。那么就要用到makeblastdb。这个命令很简单。USAGE makeblastdb [-h] [-help] [-in input_file] [-input_type type] -dbtype molecule_type [-title database_title] [-parse_seqids] [-h

SnpEff is a variant annotation and effect prediction tool. It annotates and predicts the effects of variants on genes 详细的说明请阅读：http://snpeff.sourceforge.net/SnpEff_manual.html一， 安装：首先在家目录下， 下载安装包wget h

最近在使用RNA_seq数据做些分析，结果得到了大量差异表达以及共表达的基因，如何合理展示这些基因也是一件不简单的事情。除了常见的热图(heatmap)展现形式,今天在推荐另外一种展示方式(上图C)。需要R包TCseq或者Mfuzz。我这里给出的代码是基于TCseq。 输入文件就是一个基因表达量的矩阵，如下图。 代码也很简单,见下图，也请点击阅读原文查看代码。​ library(TCseq)

根据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')

在线绘制基因表达热图昨天需要画基因表达热图(heat map), 想起一个在线的网站可以很方便的绘制，提供大家。不用R代码，这个工具绝对可以画出你想要的那样。网址是：https://software.broadinstitute.org/morpheus/

输入文件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

我的青春     最近在做一些小麦基因的表达分析，想到使用RNA-seq的数据进行生物信息学分析，并且比我做实验用的组织还要多。     下载数据之后，首先要对数据进行低质量序列和载体序列等污染序列去除，我这里结合了两个软件AdapterRemoval和bbduk2, bbduk2是bbmap中的一个子程序。AdapterRemoval --file1 input1.fastq.gz --fil

Sambamba: process your BAM data faster!　　对于很大的(>100G)的bam文件，排序时间很慢不说，往往需要1天或更多的时间，但结果还会出错。如下边的错误. 经测试Sambamba表现较好，能够节省很多时间。随着接触的数据越来越多，感觉很简单的事情也需要花很多时间。不仅仅是数据多了的问题！[bam_sort_core] merging from 3288 fil

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

Magic-BLAST is a tool for mapping large next-generation RNA or DNA sequencing runs against a whole genome or transcriptome.前面一段 时间看到NCBI上blast出了一个新系列，短序列的mapping。基本与blast使用体验一致。现在还未测试Magic-BLAST 是否适用于长

MOSAIK安装1、下载git clone https://github.com/wanpinglee/MOSAIK.git2、安装进入src文件夹，然后make #注意如果静态库没有安装可使用 yum install glibc-static 安装#或者修改includes/linux.inc文件，将"-static"删除3、使用还未探索，这是文章MOSAIK: A Hash-Based Alg

经常遇到提取问题，一旦要提取的文件很大，关键字很多，可以使用集合#!/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

使用gmap软件输出的gff3文件并不能直接用来作为gffcompare的输入文件 故此，需要进行转换，转换也简单就是讲target字段删除即可。#!/usr/bin/env python# -*- coding: utf-8 -*- __author__ = 'shengwei ma'__author_email__ = 'shengweima@icloud.com'with open('

>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

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

负号表示反向序列，并且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

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

##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

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

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

输入文件格式：UN182534 rootUN182535 rootUN182536 rootUN182537 rootUN052459 flower seed seed seed flower flower root seed flower flower flower root root flower seed flower flower seed seed leaf flower

get fasta sequence from gff3chr1A NRGenome exon 1157233 1158291 100 + . ID=TRIAE_CS42_U_TGACv1_641506_AA2096860.1.mrna1.exon1;Name=TRIAE_CS42_U_TGACv1_641506_AA2096860.1;Parent=TRIAE_CS42_U_TGACv1_6

IntroductionWhat is HISAT2?HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome, transcriptome, and exome sequencing data) against the general h

我是生物信息学专业，研究对象是小麦，小麦是异源六倍体作物，基因组多达17Gb，是水稻基因组的40倍，大约85%是重复序列。至今，小麦基因组仍然没有一个较好的参考基因组，好在，国际小麦测序组织已经公布了一个相对较好的reference，但是对应的基因组注释没有公布。我实在仍不住吐槽，很多工作放到水稻上1个月可以解决的东西，放到小麦上需要至少1年多，当然还不一定解决。我们没有自己测序的数据，加上小麦基因

前言本博文基本是翻译自文章A beginner’s guide to eukaryotic genome annotation 　　 二代测序以及最近三代单分子测序的火热，让我们获得高质量基因组越来越来容易，然而基因组注释仍然面临许多挑战。其中一个挑战就是寻找基因（gene finding），训练基因model，选择基因预测软件和注释软件，另一个挑战就是更新合并不同途径注释的基因，目前还没有完

根据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