提取指定区间的转录本序列

用户需求：给定一个区间，找到该区间内所有的基因、转录本。
获得所有的转录本序列。
解决方案：
1.获取目标区段内所有的基因、转录本

cat merged.gtf|awk '$1 ~/8/ && $4>90193841 && $5<91589697 {print $0}' >tran_8

tran_8即是给出所有的在区间内的基因、转录本。
参数讲解：
awk $1匹配第一列 以8开头的，此处是第8染色体。 && &运算符，同时满足三个条件的，返回后面指定的列。

2.从基因组提取指定的序列。
cat tran_8 |awk '$2 ~/trans/ {print $8}|wc -l'
统计总共有多少个转录本。
实现方案：

• bedtools（常用工具）
• pysam (python包）
  先用bedtools，据说速度飞快。
  参考链接

说明一点：bedtools其中一个输入文件格式要求是tab分隔符的，我的解决方案是先用awk选好文件内容，之后导入本地的excel，另存为制表符格式。否则，bedtools会一直报错，说你的输入不够3列
我太喜欢awk,sed,grep这种神奇的命令，快捷干净，无废话。
1.先从所有的文件中找出你需要的区间内的所有转录本位置
cat merged.gtf|awk '$1 ~/8/ && $4>90193841 && $5<91589697 {print $0}' >tran_8

2. 导入本地，转换为制表符格式。注意，尽量格式顺序如下

   
   image.png
   

   第4列是你要的命名的转录本名称。把第9列和第5列合并，即为第4列。第一列是chr8还是8根据你的基因组序列里面的格式决定。
   输出文件名为tran_8.txt

3. bedtools转换参数

bedtools getfasta -fi /disks/backup/chaim/maize/Zea_mays.B73_RefGen_v4.42.fa -bed tran_8.txt -split -name -fo trans_fa

最后输出文件为trans_fa
4.提取我们需要的转录本--合并后的转录本

sed -n '/trans/{N;p}' trans_fa >transcript

transcript即为我们的最终结果。

经验：学习软件用法最准确快捷的地方是软件本身的文档或guide,各种教程的用法不一定准确方便，还浪费精力。